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ABSTRACT 

We compare Tully-Fisher (TF) data for 838 galaxies within cz = 3000 kms -1 from the Mark 
III catalog to the peculiar velocity and density fields predicted from the 1.2 Jy IRAS redshift 
survey. Our goal is to test the relation between the galaxy density and velocity fields predicted by 
gravitational instability theory and linear biasing, and thereby to estimate /?/ = fi a6 /6/, where bj is 
the linear bias parameter for IRAS galaxies on a 300 kms -1 scale. Adopting the IRAS velocity and 
density fields as a prior model, we maximize the likelihood of the raw TF observables, taking into 
account the full range of selection effects and properly treating triple-valued zones in the redshift- 
distance relation. Extensive tests with realistic simulated galaxy catalogs demonstrate that the 
method produces unbiased estimates of (3j and its error. When we apply the method to the real 
data, we model the presence of a small but significant velocity quadrupole residual (~ 3.3% of Hubble 
flow), which we argue is due to density fluctuations incompletely sampled by IRAS. The method 
then yields a maximum likelihood estimate /3j = 0.49 ± 0.07 (1 a error). We discuss the constraints 
on and biasing that follow from this estimate of /3j if we assume a COBE-normalized CDM power 
spectrum. Our model also yields the one dimensional noise in the velocity field, including IRAS 
prediction errors, which we find to be 125 ± 20 kms -1 . 

We define a x 2 -like statistic, x|, that measures the coherence of residuals between the TF data 
and the IRAS model. In contrast with maximum likelihood, this statistic can identify poor fits, 
but is relatively insensitive to the best /?/. As measured by %|, the IRAS model does not fit the 
data well without accounting for the residual quadrupole; when the quadrupole is added the fit is 
acceptable for 0.3 < (3i < 0.9. We discuss this in view of the Davis, Nusser, & Willick analysis that 
questions the consistency of the TF and IRAS data. 

1. Introduction 

One of the most important tasks facing observational cosmology is determination of the density parameter 
Q. Along with the Hubble constant Hq and the cosmological constant A, the density parameter fixes the global 
structure of spacetime. One approach to the problem uses the classical cosmological tests of the geometry 
of the universe, such as the apparent magnitudes as a function of redshift of standard candles (e.g., Type la 
Supernovae, Perlmutter et aJ. 1996). While promising, this approach is sensitive to the possible evolution of 
the standard candles with redshift. Moreover, it is difficult to disentangle the effects of A and in such tests 
(Dekel, Burstein, & White 1997). Alternatively, one may carry out dynamical measurements of f2 in the local 
(z ^ 0.05) universe, in which both evolution and the geometrical effects of the cosmological constant may be 
safely neglected. 

Low-redshift tests of f2 are based on dynamical measurements of the mass of gravitating matter on some 
characteristic size scale. For example, measurements of rotation curves (Rubin 1983) or the motions of satellite 
galaxies (Zaritsky et ai. 1993) yield the masses of ordinary spirals within ~ 10-200 kpc of their centers. The 
velocity dispersions (Carlberg et al. 1996), X-ray temperatures (White et al. 1993), and gravitational lensing 
effects (Tyson & Fischer 1995; Squires et al. 1996) of rich clusters of galaxies provide mass estimates on ~ 1 Mpc 
scales. In general, these and other dynamical analyses of matter in the highly clustered regime have pointed to 
a mass density corresponding to f2 ~ 0.2 ± 0.1 (e.g., Bahcall, Lubin, & Dorman 1995). This value exceeds that 
implied by known sources of luminosity (Qi nm ^ 0.01; Peebles 1993) or inferred from primordial nucleosynthesis 
(f^baryon ~ 0.05; Turner et al. 1996), and thus points to the existence of nonbaryonic dark matter. However, 
it is well below the Einstein-de Sitter value of CI = 1 that is favored by simplicity and coincidence arguments 
(e.g., Dicke 1970). The natural expectation from the inflation scenario is that the universe is flat, 17 + = 1, 
where Q\ = A/3-f/g is the effective energy density contributed by a cosmological constant (Guth 1981; Linde 
1982; Albrecht & Steinhardt 1982). However, if Q ~ 0.2, this inflationary prediction requires f^A — 0.8, which 
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conflicts with upper limits obtained from studies of gravitational lensing (Carroll, Press, & Turner 1992; Maoz 
& Rix 1993; Kochanek 1996). 

It is possible, however, that £1 could be close to or exactly equal to unity despite evidence to the contrary 
from dynamical tests on ~ 1 Mpc scales. This could occur if the dark matter is poorly traced by dense 
concentrations of luminous matter such as galaxies and galaxy clusters. If so, dynamical tests on scales ^ 10 
Mpc are necessary to obtain an unbiased estimate of O. Such tests involve measurements of the coherent, large- 
scale peculiar velocities of galaxies. According to gravitational instability theory (cf. Eq. [I]), these motions are 
related in an ^-dependent way to the large-scale distribution of mass. If the latter, in turn, can be inferred 
from the observed distribution of galaxies on large scales, one might hope to derive an estimate of f2 that is 
free from the pitfalls of small-scale dynamical analyses. 

This program requires a comparative analysis of two types of data sets. The first consists of radial velocities 
and redshift-independent distance estimates for large samples of galaxies. The largest such compilation to date 
is the Mark III catalog (Willick et al. 1997), which contains distance estimates for ~ 3000 spiral galaxies from 
the Tully-Fisher (1977; TF) relation, and for 544 elliptical galaxies from the D n -a relation (Djorgovski &; Davis 
1987; Dressier et al. 1987). The second type of data set is a full-sky redshift survey with well-understood 
selection criteria. Several large redshift surveys exist (cf. Strauss & Willick 1995, hereafter SW, and Strauss 
1996a, for reviews); the one which most nearly meets the requirements of full-sky coverage and well-understood 
selection is the IRAS 1.2 Jy survey (Fisher et al. 1995). The basic idea behind the comparison is as follows. 
In the linear regime (mass density fluctuations 5 = 5p/po <C 1), the global relationship between the peculiar 
velocity field v(r) and the mass-density fluctuation field 5(r) is given by gravitational instability theory: 

4 7r J |r' — r| 

where /(fi) « Q os (Peebles 1980) .Q If mass density fluctuations are equal to galaxy number density fluctuations, 
at least on the scales ( ^ few Mpc) over which it is possible to define continuous density fields, then the redshift 
survey data yield a map of <5(r) (after correction for peculiar velocities; Appendix |A|). By Eq. (||), one then 
derives a predicted peculiar velocity field v(r) as a function of The TF or D n -a data provide the observed 
peculiar velocities. The best estimate of £1 is the one for which the predicted and observed peculiar velocities 
best agree. 

Two obstacles make this comparison a difficult one. The first, already alluded to, is fundamental: one 
observes galaxy number density (5 g ) rather than mass density (<5) fluctuations. A model is required for relating 
the first to the second. The simplest approximation is linear biasing, 

S g (r) = bS(r), (2) 

in which the bias parameter b is assumed to be spatially constant. Substituting Eq. (|J) in Eq. (|l|) yields 

f 3 , a g (r')(r'-r) 
V(r) = 4^7 dr |r'-r|3 ' (3) 

where (3 = f{Q)/b. Thus, under the dual assumptions of linear dynamics and linear biasing, comparisons of 
peculiar velocity and redshift survey data, by themselves, can yield the parameter (3 but not Q. One might 
hope to break the Q-b degeneracy by generalizing Eq. (Q) to the nonlinear dynamical regime (cf. Dekel 1994, 
§ 2, or Sahni & Coles 1996, for a review). However, such generalizations are difficult to implement in practice; 
furthermore, nonlinear extensions to Eq. (H) will enter to the same order as nonlinear dynamics (we discuss 



1 We measure distances r in velocity units (km s 1 ). In such a system of units, the Hubble Constant is equal to unity by definition, and does not affect the 
amplitude of predicted peculiar velocities. 
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this issue further in § 3.3.1). Thus, without a more realistic a priori model of the relative distribution of 
galaxies versus mass, it is prudent to limit the goals of the peculiar velocity-redshift survey comparison to 
testing gravitational instability theory and determining (3. One may then adduce external information on the 
value of b to place constraints on 0, itself. 

The second obstacle is the sheer technical difficulty of the problem. The redshift-independent distances 
obtained from methods such as TF are large (~ 20%; Willick et al. 1996), and are subject to potential systematic 
errors due to statistical bias effects (Dekel 1994; SW, § 6). Furthermore, we measure the galaxy density field 5 g 
in redshift space, whereas it is the real-space density that yields peculiar velocities via Eq. (JJ). The relationship 
between the two depends on the peculiar velocity field itself. Self-consistent methods, in which v p is both the 
desired end product and a necessary intermediate ingredient in the calculation, must therefore be developed 
for predicting peculiar velocities from redshift surveys (Appendix [A]). For these reasons, reliable comparisons 
of peculiar velocity and redshift survey data require extremely careful statistical analyses. 

This problem has inspired a number of independent approaches in recent years. The potent method (Dekel, 
Bertschinger, & Faber 1990; Dekel 1994; Dekel et al. 1997) was the first effort at a rigorous treatment of 
peculiar velocity data. Dekel et al. (1993) compared the potent reconstruction of the Mark II peculiar velocity 
data (Burstein 1989) to the IRAS 1.936 Jy redshift survey (Strauss et al. 1992b), finding /?/ = 1.28±g;^| at 
95% confidence.^ An improved treatment using the Mark III peculiar velocities (Willick et al. 1997) and the 
IRAS 1.2 Jy survey (Fisher et al. 1995) yields /?/ = 0.86 ± 0.15 (Sigad et al. 1997, hereafter potiras). Hudson 
et al. (1995) compared the optical redshift survey data of Hudson (1993) to the potent reconstruction based 
on a preliminary version of the Mark III catalog, finding /3 opt = 0.74 ± 0.13 (lcr errors). These results from 
potent were obtained using 1200 kms" 1 Gaussian smoothing. A distinct approach, which differs from potiras 
in the statistical biases to which it is vulnerable (SW), and which typically uses much smaller smoothing, is 
to predict galaxy peculiar velocities and thus distances from the density field, and then use these predictions 
to minimize the scatter in the TF or D n -o relations (Strauss 1989; Hudson 1994; Roth 1994; Schlegel 1995; 
Shaya, Peebles, &: Tully 1995; Davis, Nusser, & Willick 1996, hereafter DNW). This second kind of analysis 
has produced estimates of (3i in the range ~ 0.4-0.7, lower than the values obtained from potiras. We further 



clarify the distinction between the two methods in § 2.1, and discuss possible reasons for the discrepancies in 



In this paper, we present a new maximum-likelihood method for comparing TF data to the predicted peculiar 
velocity and density fields in order to estimate (3. Its chief strength is an improved treatment of nearby galaxies 
(cz < 3000 kms -1 ), and we limit the analysis to this range. The TF data we use comprise a subset of the Mark 
III catalog of Willick et al. (1997). The predicted peculiar velocities are obtained using new reconstruction 
methods (Appendix ^) from the IRAS 1.2 Jy redshift survey.^ The outline of this paper is as follows. In 
§ ||, we first review the strengths and weaknesses of existing approaches, and then describe our new method in 
detail. In § ||, we present tests of the method using mock catalogs. In § ||, we apply the method to the Mark 
III catalog and obtain an estimate of (3j. In § ||, we analyze residuals from our maximum likelihood solution 
in order to assess whether IRAS predictions give a statistically acceptable fit to the Mark III data. In § ||, 
we further discuss and summarize our principal results. This paper is the product of nearly three years work 
and contains considerable detail. We recommend that readers interested primarily in results and interpretation 
skim § 2, and then read § 3.1, 4.4, 4.5, 5.1, 5.2, and 6.4. 



2 Because the bias parameter can differ for different galaxy samples, the value of /3 can differ as well. We will use f3j for the IRAS redshift survey and /3 pt f° r an 
optical survey. Because optical galaxies are about 30% more clustered than IRAS galaxies (SW), the conversion is /3j ~ 1.3/3 pt- When speaking gcncrically about 
the velocity-density relation, we will place no subscript on /3. 

^The original IRAS 1.936 Jy survey was presented in a series of six papers (Strauss et al. 1990; Yahil ct al. 1991; Davis, Strauss &z Yahil 1991; Strauss ct al. 1992abc), 
numbered 1, 2, 3, 4, 5, and 7, respectively. The missing paper 6 was to be the comparison of the observed and predicted velocities, to be based on Chapter 3 of 
Strauss (1989). However, it has taken us until now to come up with statistically rigorous ways of doing this comparison. The long-lost IRAS Paper 6 has thus been 
incorporated into Dekel ct al. (1993), DNW, Sigad ct al. (1997), and especially this paper. 
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2. Description of the Maximum Likelihood Method 
2.1. Alternative Approaches to the Peculiar Velocity-Density Comparison 

Before presenting our method in detail, we briefly review the principal alternatives. Two approaches are 
fairly paradigmatic, and serve to illustrate the main issues and motivate our approach. These are the potent 
method of Dekel and coworkers (e.g., Dekel 1994; Dekel et al. 1997) mentioned in § [l], and the itf method of 
Nusser & Davis (Nusser & Davis 1995; DNW). 

The potent algorithm is designed to reconstruct, from sparse and noisy radial peculiar velocity estimates, 
a smooth three-dimensional peculiar velocity field and the associated mass density field. The method is based 
on the property that the smoothed velocity field of gravitating systems is the gradient of a potential. The 
divergence of Eq. (||) is 

V • v = -P5 g . (4) 

Thus, (3 is the slope of the correlation between V • v, obtained from potent, and 5 g obtained from redshift 
survey data. This is the basis of the potiras approachP] to determining fij, discussed above (§ [j]). 

potent has several advantages as a reconstruction method. It yields model-independent, three-dimensional 
velocity and density fields well-suited for comparison with theory and for visualization. It works in the space 
of TF-inferred distances, i.e., it is a Method I approach to velocity analysis (cf. SW, § 6.4.1). Unlike Method II 
approaches (see below), it does not assume that there is a unique distance corresponding to a given redshift. In 
regions where galaxies at different distances are superposed in redshift space, potent is capable of recovering 
the true velocity field. The potiras comparison between the mass and galaxy density fields is entirely local 
(Eq. ||]), whereas predicted peculiar velocities are highly nonlocal (Eq. ||). Locality ensures that biases due to 
unsampled regions are minimized. 

The liabilities of potent are closely related to its strengths. In order to construct a model-independent 
velocity field it must have redshift-independent distances as input. Such distances require properly calibrated 
TF relations. In particular, the TF distances for samples that probe different regions of the sky must be 
brought to a uniform system, which is a difficult procedure (cf. Willick et al. 1995, 1996, 1997). Errors made in 
calibrating and homogenizing the TF relations will propagate into the potent velocity field. Because potent 
works in inferred distance space, it is subject to inhomogeneous Malmquist bias (Dekel, Bertschinger, & Faber 
1990). Minimizing this bias requires significant smoothing of the input data, potent currently employs a 
Gaussian smoothing scale of 1000-1200 kms -1 (Dekel 1994; Dekel et al. 1997), making it relatively insensitive 
to dynamical effects on small scales. As a result, the current potent applications are not particularly effective 
at extracting detailed information from the velocity field in the local (cZrS 3000 kms" 1 ) universe. 

DNW take a different approach. They work with the "inverse" form of the TF relation (Dekel 1994, § 4.4; 
SW, § 6.4.4), and thus refer to their method as itf. They express peculiar velocity as a function of redshift- 
space, rather than real-space, position; in the terminology of SW, itf is thus a Method II analysis, largely 
impervious to inhomogeneous Malmquist bias. DNW expand the redshift-space peculiar velocity field in a set 
of independent basis functions, or modes, whose coefficients are solved for simultaneously with the parameters of 
a global inverse TF relation via x 2 minimization of TF residuals. The TF data are never converted into inferred 
distances and thus do not require pre-calibrated TF relations. The IRAS-predicted velocity field is expanded 
in the same set of basis functions, allowing a mode-by-mode comparison of predicted and observed peculiar 
velocities. This ensures that one is comparing quantities that have undergone the same spatial smoothing, a 
desirable characteristic of the fit. 

As with potent, the strengths of itf are connected with certain disadvantages. Because it is a Method II 
approach, multivalued or flat zones in the redshift-distance relation (see below) necessarily bias the itf analysis. 



4 Dckcl ct al. (1993) and Sigad ct al. (1997) actually use a non-linear extension to Eq. 
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It neglects the role of small-scale velocity noise, which is non-negligible for galaxies within 1000 kms . These 
features make itf, like potent, a relatively ineffective tool for probing the very local region. Last and most 
importantly, the itf method as implemented by DNW requires that the raw magnitude and velocity width data 
from several distinct data sets be carefully matched before being input to the algorithm. Any systematic errors 
incurred in matching the raw data from different parts of the sky will induce large-scale, systematic errors in 
the derived velocity field. Thus, although itf does not need input TF distances, it is vulnerable to a priori 
calibration errors just as potent is. 



2.2. VELMOD 

The approach we take in this paper, "velmod," is a maximum likelihood method designed to surmount 
several of the difficulties that face potiras and itf. velmod generalizes and improves upon the Method II 
approach to velocity analysis. Method II takes as its basic input the TF observables (apparent magnitude and 
velocity width) and redshift of a galaxy, and asks, what is the probability of observing the former, given the 
value of the latter? It then maximizes this probability over the entire data set with respect to parameters 
describing the TF relation and the velocity field. The underlying assumption of Method II is that a galaxy's 
redshift, in combination with the correct model of the velocity field, yields its true distance, which then allows 
the probability of the TF observables to be computed. This analytic approach was originally developed by 
Schechter (1980), and was later used by Aaronson et ah (1982b), Faber & Burstein (1988), Strauss (1989), Han 
& Mould (1990), Hudson (1994), Roth (1994), and Schlegel (1995), among others. 

The main problem with Method II is its assumption that a unique redshift-distance mapping is possible. 
This assumption breaks down for two reasons. First, redshift is a "noisy" realization of distance plus predicted 
peculiar velocity — both because of true velocity noise generated on very small ( £j 1 Mpc) scales, and because 
of the inaccuracy of the velocity model (even for the correct (3) due to nonlinear effects and shot noise in the 
density field. Second, even in the absence of noise, the redshift-distance relation can, in principle, be multivalued: 
more than one distance along the line of sight can correspond to a given redshift. velmod accounts for all of 
these effects statistically by replacing the unique distance of Method II with the joint probability distribution 
of redshift and distance. This distribution is constructed to allow for both noise and multivaluedness. The 



distance dependence is then integrated out (§ |2.2.l|) , yielding the correct probability distribution of the TF 
observables given redshift. 

There are two additional advantages to the velmod approach. First, it requires neither a priori calibration 
of the TF relations (as does potent) nor matching of the input data from disparate samples (as does itf). 
An individual TF calibration for each independent sample occurs naturally as part of the analysis. Second, it 
does not require smoothing of the input TF data, and thus allows as high-resolution an analysis as the data 
intrinsically permit. This second feature, along with its allowance for velocity noise and triple-valued zones, 
makes velmod well-suited for probing the local (cz^,3000 kms -1 ) velocity field. An analysis of local data is 
desirable because random and systematic errors in both the IRAS and TF data are less important nearby than 
far away. 



2.2.1. Mathematical Details 

We now describe the method in detail. We assume that the relevant distance indicator is the TF relation; 
with minor changes the formalism could be adapted to comparable distance indicators such as D n -a. We use 
the terminology of Willick (1994) and Willick et aJ. (1995): briefly, we denote by m and 77 = logu ro t — 2.5 a 
galaxy's corrected apparent magnitude and velocity width parameter, respectively; by cz its Local Group frame 
radial velocity ("redshift") in kms -1 ; and by r its true distance in kms -1 . We define the distance modulus as 
H = 5 log r, and absolute magnitudes as M = m — /x. We write the forward and inverse TF relations as linear 
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expressions, M(rf) = A — brj and rp(M) = — e(M — D), and denote their rms scatters cttf and respectively. 

We seek an exact expression for the probability that a galaxy at redshift cz possesses TF observables (m, 77) 
given a model of the peculiar velocity and density fields.^ We first consider the joint probability distribution of 
the TF observables, redshift, and an unobservable quantity, the true distance r. Later, we will integrate over r 
to obtain the probability distribution of the observables. We may write 

P(m, 77, cz, r) = P(m, rj\r) x P(cz\r) x P{r) . (5) 

The splitting into conditional probabilities reflects the fact that the TF observables and the redshift couple 
with one another only via their individual dependences on the true distance r. 

The first of the three terms on the right hand side of Eq. (||) depends on the luminosity function, the sample 
selection function, and the TF relation. We can express it in one of two ways, depending on whether we are 
using the forward or inverse form of the TF relation: 

1. Forward relation: 

p/ I \ a>{ \c( ^ 1 ( [m-(M(r/)+/i(r))] 2 \ 

P{m,ri\r) ck (p{r])b{rn,rj,r) exp — s (o) 

cttf \ 2cr TF J 

2. Inverse relation: 

P(m, r]\r) oc $(m — fj,(r))S(m, 77, r) — exp — 5 — , (7) 

a v V 2a v J 

where <j)(rj) and $(M) are the (closely related) velocity width distribution function and luminosity function, 
S(m, 77, r) is the sample selection function, and we have assumed Gaussian scatter of the TF relation (an 
assumption validated by Willick et aJ. 1997). Detailed derivations of these expressions are given by Willick 
(1994).[] In Eqs. (||) and (|7j) we have written only proportionalities, as the normalization is straightforward and 
will occur at a later point in any case. 

The third term on the right hand side of Eq. (||) is simply the a priori probability of observing an object at 
distance r, 

P(r) oc r 2 n(r) , (8) 

where n(r) oc 1 + S g (r) is the number density of the species of galaxies that makes up the sample. The second 
term on the right hand side of Eq. (|5|), P(cz\r), is the one which couples the TF observables to the velocity 
field model. We assume that, for the correct IRAS velocity field reconstruction (i.e., for the correct value of 0i 
and other velocity field parameters to be described below), the redshift is normally distributed about the value 
predicted from the velocity model: 

, x 1 ( \cz- (r + u(r))} 2 \ , . 

P(czr) = ^=^exp -i \ 1 ,n , 9 

where u(r) = r • [v(r) — v(0)] is the radial component of the predicted peculiar velocity field in the Local Group 
frame (cf. Eq. |A1|) . We treat the velocity noise a v as a free parameter in our analysis; we discuss its origin 



J Thc dependence of all quantities on the line of sight direction will remain implicit. 



^Willick (1994) assumed that the selection function depended only on the TF observables. Here, we acknowledge the possibility of an explicit distance dependence; 
the origin of such a dependence was discussed by SW, § 6.5.3. 
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in detail in § 3.2. Although a v must be position or density dependent at some level, we treat it as spatially 
constant in this paper, except in the Virgo cluster (§ |4.3| ). 

Substituting Eqs. @ or |7|), (^|), and (|9|) into Eq. (|5|) yields the joint probability distribution P(m,rj,cz,r). 
To obtain the joint probability distribution of the observable quantities, one integrates over the (unobserved) 
line-of-sight distance, i.e., 

P(m,7],cz) = / P(m,7],cz,r) dr . (10) 
Jo 

In practice, it is not optimal to base a likelihood analysis on the joint distribution P(m, rj, cz) because of its 
sensitivity to terms, such as the luminosity function, the sample selection function, and the density field, that 
are not critical for our purposes. Instead, the desired probability distributions are the conditional ones: 

1. Forward TF relation: 

P(m, rj, cz) 



P(m\rj, cz) 



J^o P{m, rj, cz) dm 



J^°dr r 2 n(r) P(cz\r) S(m, rj, r) exp ( _ i m ( M (v)+n( r ))] 



Jo^dr r 2 n(r) P(cz\r) 'dm S(m, 77, r) exp ( _ i m ( M (^)+M r ))] 

\ Za TF 

2. Inverse TF relation: 

d/ 1 \ P(m,r],cz) 
P{r]\m,cz) — 



JZo p ( m » V, cz) dr} 

roo 1 9 / \ x.t 1 w t-./ 1 \ ni \ ( \v~ V°{ m ~ MM)!' 

J dr r n[r) <P(m — /U(rj) P(cz\r) b(m, rj, r) exp I — 1 ^ L - 



J °°dr r 2 n(r) <I>(m — /u(r)) P(cz\r) J^dn S(m, m r) exp ( — ^ v 



(12) 



where P(cz\r) is given by Eq. Although neither of these expressions is independent of the density field 
n(r) or the selection function S, their appearance in both the numerator and denominator much reduces their 
sensitivity to them. A similar statement holds for the luminosity function in Eq. (|l2|). The velocity width 
distribution function <j) has, however, dropped out entirely from the forward relation probability. We discuss 
these points further in § 2.2.2| . 



Equations (0) and (12) are the conditional probabilities whose products over all galaxies in the sample we 



wish to maximize. In practice, we do so by minimizing the quantities 

C iorw = -2^2 In P(mi\r]i,czi) (13) 



or 



Anv = - 2 ^2 lnP iVi\ m i,CZi) , (14) 

i 

where the index i runs over all objects in the TF sample. We have assumed that the probabilities for each 
galaxy are independent; we validate this assumption a posteriori (cf. § |5.2j ). 



S 



2.2.2. Further discussion of the velmod likelihood 

The physical meaning of the velmod likelihood expressions is clarified by considering them in a suitable 
limit. If we take a v to be "small," in a sense to be made precise below, the integrals in Eqs. (|ll]) and ( |T2] ) may 
be approximated using standard techniques. If in addition we neglect sample selection (S = 1) and density 
variations (n(r) = constant), and assume that the redshift-distance relation is single- valued, we find for the 
forward relation: 



P(m\r], cz) 



1 



27T£T f 



exp 



1 

2a! 



m 



M(rj) + b\ogw + 



10 
hTlO" 



A 



(15) 



where w is the solution to the equation cz = w + u(w), i.e., it is the distance inferred from the redshift and 
peculiar velocity model; A„ = a v /[w(l + u')], where u' = (du/dr) r=w , is the effective logarithmic velocity 
dispersion; and 



<7tF + 



In 10 



A 



1/2 



(16) 



is the effective TF scatter, including the contribution due to a v . An analogous result holds for the inverse 
relation. The criterion A 2 , <C 1, which quantifies the statement that o~ v is "small," must be satisfied to derive 
Eq. (H). 

Eq. (15) shows that the probability distribution P(m\rj, cz) preserves the Gaussian character of the real- 
space TF probability distribution P(m\rj,r) in this limit. However, the expected value of m is shifted from the 
"naive" value M(rj) + h\ogw by an amount ~ 4.3A 2 . This shift is in fact nothing more than the homogeneous 
Malmquist bias due to small-scale velocity noise; it differs in detail from the usual Malmquist expression (i.e., 
that which affects a Method I analysis) because it arises from the Gaussian (rather than log-normal) probability 
distribution, Eq. @. Furthermore, the effective scatter a e is larger than cttf, because the velocity dispersion 
introduces additional distance error and thus magnitude scatter. The effects associated with velocity noise 
diminish with distance (A^ oc r -1 ), however; the velocity Malmquist effect vanishes in the limit of large 
distances, in contrast with the distance-independent Malmquist effect for Method I, and the effective scatter 
approaches the TF scatter. At large enough distance the velmod likelihood approaches a simple Gaussian TF 
distribution with expected apparent magnitude M(rj) + 51ogw;, and velmod reduces to standard Method II. 

Indeed, Eq. ( |l5| ) enables us to define the regime in which velmod represents a significant modification of 
Method II. The distance ru at which the velocity noise effects become unimportant is determined by m 3> 
<t„/Atf(1 + u'), where Atf = hil0<7TF/5 is the fractional distance error due to the TF scatter (Atf — 0.2 
for the samples used here). For a v = 125 kms -1 , the value we find for the real data (§ 4.5), this shows that in 
the unperturbed Hubble flow, where u' = 0, velocity noise effects become unimportant beyond ~ 1500 kms -1 . 
However, at about this distance, in many directions, the Local Supercluster significantly retards the Hubble flow, 
u' ~ —0.5, so that the effective a v is about twice its nominal value. Thus, velmod in fact differs substantially 
from Method II to roughly twice the Virgo distance. This fact guided our decision to apply velmod only out 
to 3000 kms" 1 (cf. §@). 

Eq. (|l5|) also demonstrates that maximizing likelihood (minimizing £f or w) is not equivalent to \ 2 minimiza- 
tion, even under the adopted assumptions of constant density and negligible selection effects, because of the 
factor o~~ 1 in front of the exponential factor. This factor couples the velocity model (i.e., the values of w and 
u'(uu)) to the velocity noise. In particular, maximizing the velmod likelihood is not equivalent to minimizing 
TF scatter (cf. § |4.5| ), except in the limit that a v is set to zero. 

The assumptions required for deriving Eq. ( |T5| ) remind us that there are two other factors which distinguish 
velmod from standard Method II. First, for realistic samples one cannot assume that 5 = 1. The presence of the 
selection function in Eqs. (||) and (^) is essential for evaluating true likelihoods, and we have fully incorporated 
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these effects into our analysis.^ Second, the gaiaxy density n(r) is not effectively constant along most lines of 
sight. Thus, velmod, like Method I but unlike Method II, requires that n(r) be modeled. We do so here by 
using the IRAS density field itself, which is a good approximation to the number density of the spiral galaxies 
in the TF samples. The density field has a non- negligible effect on the velmod likelihood whenever it changes 
rapidly on the scale of the effective velocity dispersion a v /(l + u'). 

The most significant differences between velmod and Method II thus occur in regions where v! — ► —1 (flat or 
triple- valued zones), or when the density varies particularly sharply. In practice, both these effects occur in the 
vicinity of large density enhancements such as the Virgo cluster. We illustrate this in Figure 1, which shows the 
redshift-distance relation, and the corresponding value of P(r\cz) oc P(cz\r)P{r) in the vicinity of triple- valued 
zones. When looking at these panels, keep in mind that the velmod likelihood is given by multiplying P{r\cz) 
and the TF probability factor P(m\rj,r) and integrating over the entire line of sight. Panels (a) and (b) depict 
the situation near the core of a strong cluster, and panels (c) and (d) farther from the center. In each case, the 
cloud of points represents the velocity noise, here taken to be a v = 150 kms -1 . In panel (a), the redshift of 
1200 kms -1 crosses the redshift-distance diagram at three distinct distances. The quantity P(r\cz) shows three 
distinct peaks. The highest redshift one is the strongest because of the r 2 weighting in Eq. (|J). In panel (b), the 
redshift of 1700 kms -1 is such that the object just misses being triple- valued; however, the finite scatter in the 
redshift-distance diagram means that there is still appreciable probability that the galaxy be associated with 
the near-crossing at cz ~ 900 kms -1 . In panel (c), the redshift-distance diagram goes nearly flat for almost 
600 kms -1 ; a redshift that comes close to that flat zone has a probability distribution that is quite extended. 
Finally, panel (d) shows a galaxy whose redshift crosses the redshift-distance diagram in a region in which it is 
quite linear, and the probability distribution has a single narrow peak without extensive tails. 

Two final details deserve brief mention. First, the integrals over m and r\ that appear in the denominators 



of Eqs. ( |ll| ) and (12) may be done analytically for the case of "one-catalog selection" studied by Willick (1994, 
§ 4.1), which indeed applies for the samples used in this paper (Willick et al. 1996). The numerical integrations 
required to evaluate Eqs. (|il|) and ( |l2|) are thus one-dimensional only. Second, as noted above, the velocity 
width distribution function ^(rj) drops out of Eq. ([[l]), but the luminosity function $(M) does not drop out 
of Eq. (12). Thus, inverse velmod requires that we model the luminosity function of TF galaxies. This is 
an annoyance at best, and could introduce biases, if we model it incorrectly, at worst. We have thus chosen 
to implement only forward velmod in this paper. On the other hand, inverse velmod enjoys the virtue that 
inverse Method II approaches do generally: to the degree that the selection function S is independent of r] and 
r, it drops out of Eq. (|i~2]). In a future paper, we will apply the small-cit, approximation to velmod for more 
extensive samples to larger distances. For that analysis the inverse approach will be used as well. 



2.2.3. Implementation of velmod 

The probability distribution P(m\n,cz) (Eq. [ll]) is dependent on a number of free parameters, most im- 
portantly (3j. However, because /?/ enters at an earlier stage — in the reconstruction of the underlying den- 
sity and velocity fields from IRAS (Appendix A) — it is on a different footing from other parameters. Thus, 
rather than treating f3j as a continuous free parameter, velmod is run sequentially for the ten discrete values 
01 = 0.1, . . . , 1.0 for the real data, and for the nine discrete values /?/ = 0.6, . . . , 1.4 for the mock catalog 
data (§ [|).Q For each 0j, probability is maximized (£f orw is minimized) with respect to the remaining free 
parameters. These parameters are: 

1. The TF parameters A, b, and <ttf for each sample in the analysis. Here we limit the analysis to the 



^Selection effects are not specific to VELMOD per sc, however. They can and should be modeled in any Method II-like analysis. In particular, they do not 
vanish in the A. v ► limit. 

^The choice of these values of j3j was based on the need to bracket the "true" value: 1.0 in the mock catalogs, and, as it turns out, ^0.5 for the real data. 
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Fig. 1. The effects of triple-valued or flat zones. The S-shaped curves show the relation between redshift and distance along two lines of sight to a 
cluster, (a) A galaxy with a redshift of 1200 kins - 1 can lie at three distinct distances. When the small-scale noise inherent in any velocity field model, 
as indicated by the scattered points, is taken into account, the quantity P(r\cz), shown as the three-pcaked curve at the bottom, gets smoothed out. 
(b) A galaxy at 1700 kms -1 along the same line of sight intersects the rcdshift-distance curve at only one point, but comes close enough to it elsewhere 
to give a second peak to the P(r\cz) curve, (c) Further from the cluster, the redshift-distancc curve becomes flat, giving a broad peak to P(r\cz). (d) 
At a redshift sufficiently far from a triple- valued zone, P(r\cz) has only one narrow peak. 



Mathewson et al. (1992; MAT) and Aaronson et al (1982a; A82) samples, as we discuss in § 4.1. Thus, 
there are a total of 6 TF parameters that are varied. Note that the TF scatters are not simply calculated 
a posteriori. The statistic £f orw depends on their values and they are varied to minimize it. 

The small-scale velocity dispersion o~ v . The quantities a v and cjtf can trade off to a certain extent (cf. 



Eq. 15). However, their relative importance depends on distance. Sufficiently nearby (^1000 kms" 1 ), 
a v is as large or a larger source of error than the TF scatter itself. Thus it is determined in this local 
region. Beyond ~ 2000 kms -1 , the TF scatter dominates the error, and it is determined at these distances. 
Because the samples populate a range of distances, the two can be determined separately, with relatively 
little covariance. 

3. We also allow for a Local Group random velocity vector wlg- The IRAS peculiar velocity predictions are 



given in the Local Group frame (Eq. Al). That is, the computed Local Group peculiar velocity vector 
has been subtracted from all other peculiar velocities. However, just as we expect all external galaxies to 
have a noisy as well as a systematic component to their peculiar velocity, so we must expect the Local 
Group to have one as well, especially considering the uncertainties in the conversion from heliocentric to 
Local Group frame. We allow for this by writing u(r) = uiRAs( r ) — wlg ■ r, where uiRAs(f) is given as 
described in Appendix [A], and the three Cartesian components of wlg are varied in each velmod run at a 
given /3j . We note briefly that this procedure is self-consistent only as long as |wlg| is a t most comparable 
to o~ v . In practice, we will find that for /3j near its best value, the amplitude of wlg is trivially small. 
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4. Finally, we allow for the existence of a quadrupole velocity component that is not included in the IRAS 
velocity field. The justification for such a velocity component will be discussed in § iA and Appendix B. 
The quadrupole is specified by five independent parameters, although we will not take them as free in the 
final analysis (we discuss this further in § |4.4Q . 

Thus there are 3x2 + 1 + 3 = 10 free parameters that are varied for any given value of /3j, when the 
quadrupole is held fixed. Thus, for any value of 0i, we give the data the fairest chance it possibly has to fit the 
IRAS model. In particular, the TF relations for the two separate samples used are not "precalibrated" in any 
way. This ensures that TF calibration in no way prejudices the value of /3j we derive. 



3. Tests With Simulated Galaxy Catalogs 

In this section, we test the velmod method on simulated data sets. Kolatt et al. (1996) have produced 
simulated catalogs that mimic the properties of both the IRAS redshift survey and the Mark III samples. We 
briefly review the salient points here. 

The mass density distribution of the simulated universe is based on the distribution of IRAS galaxies in the 
real universe. This was achieved by, first, taking the present redshift distribution of IRAS galaxies and solving 
for a 500 kms -1 smoothed real-space distribution via an iterative procedure that applies nonlinear corrections 
and a power-preserving filter (Sigad et al. 1997). The smoothed, filtered IRAS density field was then "taken 
back in time" using the Zekdovich-Bernoulli algorithm of Nusser & Dekel (1992) to obtain the linear initial 
density field. The method of constrained realization (Hoffman & Ribak 1991; Ganon & Hoffman 1993) was 
used to restore small-scale power down to galactic scales. The resulting initial conditions were then evolved 
forward as an = 1 iV-body simulation using the PM code of Gelb & Bertschinger (1994). The present-day 
density field resulting from this procedure is displayed in Figure 6 of Kolatt et al. (1996). 

We generated a suite of 20 mock Mark III and mock IRAS catalogs from this simulated universe]^] Each 
mock Mark III TF sample was constructed to mimic the distribution on the sky and in redshift space of the 
corresponding real sample, and the TF relations and scatters of the mock samples were chosen to be similar 
to the observed ones. The mock TF samples were subject to selection criteria similar to those imposed on 
the real samples. The mock IRAS redshift catalogs were generated so as to resemble the actual IRAS 1.2 Jy 
redshift survey. They have the true IRAS selection and luminosity functions applied, and lack data in the IRAS 
excluded zones (cf. Strauss et al. 1990). These data were then put through exactly the same code to derive 
peculiar velocity and density fields as is used for the real data (Appendix [A]). To simplify interpretation of the 
mock catalog tests, the mock IRAS galaxies were generated with probability proportional to the mass density 
itself. Thus, the mock IRAS galaxies are unbiased relative to the mass; i.e., for the mock catalogs bi = 1, and 
therefore the true value of Pi for the simulated data is unity. 



3.1. Accuracy of the (3 Determination 

The IRAS velocity field reconstructions may be produced using a variety of smoothing scales, and we have 
used 300 and 500 kms -1 Gaussian smoothing. We found, however, that at 500 kms" 1 smoothing velmod 
returned a mean Pi biased high by ~ 20%; the predicted peculiar velocities were too small, and a too-large Pi 
was needed to compensate. Our discussion from this point on will refer to 300 kms -1 smoothing, which, as we 
now describe, we found to yield correct peculiar velocities and an unbiased estimate of Pi. 



The 20 catalogs (of both types) are different statistical realizations of the same simulation. As a result, our simulations fully probe the effects pf statistical 
variance (due to distance indicator scatter, spatial inhomogcncitics. etc.) but do not include those of cosmic variance. However, as we shall argue in § |6l we expect 
that cosmic variance will have minimal effect on our /3-dctcrmination. 
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ft ft 

Fig. Z. Plots of the likelihood statistic, £f or w: versus (3i for VELMOD runs using four of the mock catalogs. (The true value of f3j for the mock 

catalogs is unity, as discussed in the text.) Also indicated on the plots arc /3 m in, the maximum likelihood values for /3j, and the average of its two 
one-sided errors 8/3±. The solid lines drawn through the points arc the cubic fits used to determine /3 m i n and 8{3±. 

velmod was run on the 20 mock catalogs, and likelihood (£forw) versus /3j curves were generated for each. 
As with the real data (§ ||), we used only the A82 and MAT TF samples; we limited the analysis to cz < 
3500 kms" 1 ™ The curves were fitted with a cubic equation of the form 

£ = £ + q(Pl- Pminf +p(Pl~ /?min) 3 (17) 

to determine (3 m { n , the value of (3j for which Cf orw is minimized. This is the maximum likelihood value of /?/. 
Four representative £f or w versus /?/ plots are shown in Figure 2, along with the cubic fits. We estimate the 1 a 
errors 5(3± in our maximum likelihood estimate by noting the values (3 ± 5(3± at which C = Cq + 1. Given the 
presence of the cubic term in Eq. (0), this is not necessarily rigorous, but we can test our errors by defining 
the x 2 -hke statistic 

X 2 = E " !) 2 /^± , (18) 

where 5(3+ was used if (3 m \ a < 1, and 5(3- was used if (3 m \ n > 1. For the 20 mock catalogs, it was found that 
X 2 = 21.2. Thus, our tests were consistent with the statement that the error estimates obtained from the change 
in the likelihood statistic near its minimum are true 1 a error estimates. Although we formally derive two-sided 
error bars, the upper and lower errors differ little, and when we discuss the real data (§ we will give only 
the average of the two. The weighted mean value of (3 m \ n over the mock catalogs was 0.984, with an error in 
the mean of ~ O.OS/y 7 ^) = 0.017. Thus, the mean (3 m - m is within ~ 1 a from the true answer. We conclude 



ll ^The real data analysis extended only to 3000 km s 1 , but because there are fewer nearby TF galaxies in the mock catalogs, we extended the mock analysis to a 
slightly larger distance. 
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that there is no statistically significant bias in the velmod estimate of /3j. The results of this and other tests 
we carried out using the mock catalogs are summarized in Table jy. 

3.2. Accuracy of the Determination of a v and wlg 

The mock catalogs also enable us to determine the reliability of the small-scale velocity dispersion a v derived 
from velmod. This quantity may be viewed as the quadrature sum of true velocity noise (<r") and IRAS velocity 
prediction errors (<7^) resulting from shot-noise and imperfectly modeled nonlinear ities. (For the real data, there 
is an additional contribution from redshift measurement errors, which are zero in the mock catalog.) We can 
measure both cr™ and a\, directly from the mock catalogs. To measure velocity noise, we determined a Ul the 
rms value of pair velocity differences oz(rj) — cz(rj) of mock catalog TF galaxies within 3500 kms -1 outside of 
the mock Virgo core, for |rj — r^| < r max . We found a u to be insensitive to the precise value of r max provided 
it was £ 150 kms -1 , implying that we are not including the gradient of the true velocity field on these scales. 
Taking 

r msLx — 150 kms \ we found <r u — 71 kms corresponding to <r™ — cr u /y2 — 50 kms This value is 
so small because the PM code does not properly model particle-particle interactions on small scales. 

We measured the IRAS prediction errors al as follows. For each mock TF particle (again, within 3500 kms -1 
and outside the mock Virgo core), we computed an IRAS predicted redshift cz\ = ri + u(ri) + fri — wlg ■ n i> 
where r, was the true distance of the object, u(rj) was the IRAS-predicted radial peculiar velocity in the Local 
Group frame (for f3j = 1), / was a zero-point error in the IRAS model (cf. § |3.3|) , and wlg was the mock 
Local Group peculiar velocity, which (just as in the real data) is not known precisely and was also treated as 
a free parameter. We then minimized the mean squared difference between cz\ and the actual redshifts czi 
over the entire TF sample with respect to / and wlg- The rms value of (czi — czf) at the minimum was then 
our estimate of the quadrature sum of IRAS prediction error and true velocity noise, which we found to be 
98 ±2 kms" 1 after averaging over the 20 mock catalogs. Subtracting off the small value of <t" found above, we 
obtain al ~ 84 kms -1 . This surprisingly small value is indicative of the high accuracy of the IRAS predictions 
for nearby galaxies not in high density environments. 



The value a v = 98 kms -1 is somewhat smaller than the real universe value of a v = 125 kms -1 (§ 4.5). 
Because we wanted the mock catalogs to reflect the errors in the real data, we added artificial velocity noise 
of 110 kms -1 to the redshift of each mock TF galaxy before applying the velmod algorithm, increasing a v to 
147 km s 1 -E3 ^he mean value of a v from the velmod runs on the 20 mock catalogs was (a v ) = 148.7 ±4.6 km s _1 , 
in excellent agreement with the expected value. We conclude that velmod produces an unbiased estimate of 
the a v , just as it does of (3j. The rms error in the determination of a v from a single realization is ~ 20 kms -1 . 

The calculation in which we minimized (czi—czf ) 2 also yielded estimates of the Cartesian components of Local 
Group random velocity vector wlg- Their mean values over 20 mock catalogs is given in Table [l], together with 
the corresponding mean values returned from velmod over the 20 mock catalog runs. The two are in impressive 
agreement. These values reflect an offset between the CMB to LG transformation assigned to the simulation 
and the average value of wlg assigned by the mock IRAS reconstruction for fii = 1. We conclude that velmod 
properly measures the Cartesian components of wlg to within ~ 50 kms -1 accuracy per mock catalog. 

3.3. The TF Parameters Obtained from velmod 

The mock catalogs also enable us to test the accuracy of the TF parameters determined from the likelihood 
maximization procedure. The comparison of input and output values is given in Table |l[ The results for the TF 
slope and scatter are consistent with the statement that velmod returns unbiased values of these TF parameters 



In retrospect, we added more noise than was necessary, but at the time we had a higher estimate of the real universe a v 
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Fig. O. Some of the parameters obtained from running VELMOD on a single mock catalog. The left hand panel shows the likelihood statistic 

along with the cubic fit used to determine its minimum. The right hand panels show the amplitude of the Local Group random velocity vector, the 
velocity noise a v , and the TF scatters for the mock A82 and mock MAT samples. Note that the Local Group velocity vector has its minimum amplitude 
for {3j ~ 1. Note also that the TF scatters do not track the likelihood curve, primarily because the velocity noise cr v also measures the inaccuracy of 
the fit. This demonstrates that minimizing TF scatter is not equivalent to maximizing likelihood, as it is for standard Method II. 

for each of the two samples. The fact that the TF scatters and a v are unbiased means that velmod correctly 
measures the overall variance in peculiar velocity predictions. 

The TF zero points returned by velmod are systematically in error by 2-3 standard deviations. This can be 
traced to a bias in the IRAS-predicted peculiar velocities; the mean value of the quantity / in § 3^ over 20 
realizations was 0.018 ± 0.007. This bias makes the IRAS-predicted distances cIiras — cz — uiras too large by 
a factor of 1.018, or ~ 0.04 mag. To bring the TF and IRAS distances into agreement, the TF zero points must 
decrease by this amount, which in fact they do (cf. Table |l]). Thus, velmod determines the TF zero points in 
such a way as to compensate for a small systematic error in the IRAS predictions. We expect such an error to 
be present in the real data as well, but it will be completely absorbed into the TF zero points, and our derived 
value of (3i will be unaffected. 



3.4. Properties of the velmod Likelihood 

The mock catalogs may also be used to illustrate some important features of the velmod analysis. An example 
of these is shown in Figure 3. The left hand panel shows Cf orw versus (3j for one of the 20 catalogs. The right 
hand panels show how three other quantities vary with /3j in the same velmod run: the amplitude of the LG 
random velocity wlg (top panel), the velocity noise a v , and the TF scatter cttf for each of the two mock TF 
samples (A82 and MAT) considered (bottom panel). Note first that the amplitude of the LG velocity vector 
is minimized near the true value of This was generally seen in the mock catalogs; it reflects the fact that 
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the fits at the wrong values of /3j try to compensate for wrong peculiar velocity predictions with Local Group 
motion. If wlg were held fixed at its maximum likelihood value, or set equal to zero, the £f OTW versus f3i curves 
would have sharper minima and the /3-uncertainty would be reduced (cf. § 4.5). Unfortunately, we cannot 
do this for the real universe because we do not know wlg a priori. Nevertheless, there we will find similar 
behavior; wlg has a minimum near the best-fit value of /?/ for the real universe. 

The figure also shows that o~ v is a weak function of /?/, but goes to a minimum at Pi ~ 1. Its value at the 
minimum for this realization is 127 kms -1 , within 1 a of the correct value (Table ||). The <ttf are also weak 
functions of but are in good agreement with the input values near the maximum likelihood values of (3j. Most 
importantly, the figure demonstrates that maximizing likelihood does not necessarily correspond to minimizing 
TF scatter, as we argued in § p. 2. 2 . The TF scatters in fact decrease monotonically with increasing /?/. As they 
do, a v increases to compensate. However, one cannot simply minimize o~tf, or even a simple combination of 
cjtf and o~ v , to obtain an unbiased (3j. One must instead maximize likelihood as defined in § [2.2.1 . 



4. Application to the Mark III Catalog Data 
4.1. Sample Selection 

To apply velmod to the real TF data, we needed first to identify a suitable subsample of the Mark III catalog. 



As discussed in § 2.2.2 , we elected to restrict the TF sample to czlg < 3000 kms , where here and throughout, 
we correct heliocentric redshifts to the Local Group frame following Yahil, Tammann, & Sandage (1977). We 
thus use the Aaronson et al. (1982a; A82) and Mathewson et al. (1992; MAT) TF samples, which are rich in 
local galaxies, restricted to this redshift interval. The two cluster samples in the Mark III catalog, HMCL and 
W91CL, contain only clusters at greater redshifts and thus are not used here. Finally, only a small fraction 
of the W91PP and CF samples is found at cz LG < 3000 kms" 1 (~ 2% of W91PP, ~ 15% of CF). This small 
number of additional galaxies was not worth the additional 6 free parameters that would be required for the 
likelihood maximization procedure (§ |2|). 

We made several further cuts on the data, as follows: 

1. An RC3 .B-magnitude limit of = 14.0 mag was adopted for A82. As discussed by Willick et al. 
(1996), A82 galaxies within 3000 kms -1 , and subject to this magnitude limit, are well described by the 
"one-catalog" selection function of Willick (1994) that enters into the likelihood equations (§ |2|). 

2. An ESO .B-band diameter limit of g?eso = arcminute was adopted for MAT. As discussed by Willick 
et al. (1996), this allows the MAT subsample to be described by the one-catalog selection function of 
Willick (1994). 

3. Only galaxies with axial ratios log(a/6) > 0.1 were included. This cut, corresponding to an inclination 
limit of % > 38° (Willick et al. 1997), reduces TF scatter due to velocity width errors. 

4. Galaxies with n < —0.45 (rotation velocities less than about 55 kms -1 ) were excluded. In practice, this 
criterion applied only to the MAT sample, which contains numerous very low-linewidth galaxies. The 
need for excluding such objects was discussed by Willick et al. (1996). 

5. Two objects within the Local Group, defined as having raw forward TF distances < 100 kms -1 , were 
excluded. No lower bound was placed on the redshifts of sample objects, however. 

This left a sample of 856 A82 and MAT galaxies. As discussed by Willick et al. (1996, 1997), real samples 
exhibit a mainly Gaussian distribution of TF residuals, but with an admixture of a few percent of non-Gaussian 
outliers. We excluded eighteen additional galaxies (4 in A82, 14 in MAT), or ~2% of our sample, because of 
their extremely large residuals from the TF relation. Finally, then, 838 galaxies, 300 in A82 and 538 in MAT, 
were used in the velmod analysis. Of these, 53 are objects found in both samples (though with different raw 
data), and thus are used twice in the analysis. 
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4.2. Velocity- Width Dependence of the TF Scatter 

It has been noted by a number of authors (Federspiel, Sandage, & Tammann 1994; Willick et aJ. 1997; 
Giovanelli et al. 1997) that cttf exhibits a velocity- width (or, equivalently, a luminosity) dependence: luminous, 
rapidly rotating galaxies have smaller TF scatter than faint, slowly rotating ones. Willick et al. (1997) showed 
that this effect could be parameterized by o~tf(v) = (J Q~9 r ]^ with different values of o~q (the scatter for a typical, 
r] = 0, galaxy) and g for each sample. For the MAT sample, they found g = 0.33, while for A82 they found 
g = 0.14. In the velmod analysis, we treated o~q as a free parameter for both samples, but fixed the values of g to 
the Willick et al. (1997) values. For the remainder of this paper, when we refer to cttf we are actually referring 
to the do for the respective samples. We note that a significant likelihood increase was achieved by adopting 
this variable TF scatter, but that the derived value of Pi was essentially unchanged. The mock catalogs were 
generated and analyzed with g = 0. 

4.3. Treatment of Virgo 

To simplify the analysis, we have taken the small-scale velocity noise, o~ v , to be independent of position. 
Clearly, this assumption must fail in the immediate vicinity of a rich cluster. The Virgo cluster is the only 
rich cluster within 3000 kms" 1 . Thus, we must artificially "cool down" the galaxies near Virgo. We do so as 
follows: if a galaxy lies within 10° of the Virgo core (taken to be I = 283.78° , b = 74.49° ) on the sky, within 
1500 kms" 1 of its mean Local Group redshift (taken to be 1035.1 kms , following Huchra 1985), and has a 
raw TF distance from Willick et al. (1997) between 800 and 2100 kms" 1 , its Local Group redshift is set to 
the mean Virgo value. Twenty objects used in the velmod analysis meet these criteria. We similarly collapsed 
mock Mark III objects associated with the Virgo cluster in the mock catalogs. 

4.4. Implementation of a Quadrupole Flow 

In the discussion of velmod in § ^, it was assumed that the IRAS-predicted velocity field, for the correct 
value of f3j, is as good a model as can be obtained. However, there can be additional contributions to the local 
flow field from structures beyond the volume surveyed (R < 12,800 kms" 1 ), as well as from shot noise- and 
Wiener-filter- induced differences between the true and derived density fields beyond 3000 kms" 1 but within 
the IRAS volume (cf. Appendix B). 

Fortunately, the nature of this contribution is such that we can straightforwardly model its general form, 
and thus treat it as a quasi- free parameter (see below) in the velmod fit. Let us write the error in the IRAS- 
predicted velocity field due to incompletely sampled fluctuations as v err (r). Because the total peculiar velocity 
field, v + v err must satisfy Eq. (Q), and because v does so by construction (Eq. ||), it follows that v err must 
have zero divergence. Moreover, if we suppose that v err corresponds to the growing mode of the linear peculiar 
velocity field, it must have zero curl as well. These properties will be satisfied if v err is given by the gradient of 
a velocity potential eft that satisfies Laplace's equation. Such a potential may be expanded in a multipole series, 
each term of which vanishes at the origin (where, by construction, v crr must itself vanish). 

The leading term in the resulting expansion of v err is a monopole, Ve2 (r) = Ar, or Hubble flow-like term. 
However, such a term is degenerate with the zero point of the TF relation (§ [3^), and is thus undetectable. The 
next term in the expansion is a dipole, virr = B, or bulk flow independent of position. Like the monopole term, 
however, the dipole term is undetectable, because we work in the frame of the Local Group. Whatever bulk 
flow is generated by distant density fluctuations is shared by the Local Group as well. The leading term in the 
expansion of v err (r) to which our method is sensitive is therefore a quadrupole term. Such a term represents the 
tidal field of mass density fluctuations not traced by the IRAS galaxies. We may write the quadrupole velocity 
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the quadrupole flow increases linearly with distance at a given position on the sky. The maximum amplitude of the quadrupole at this distance is 147 
kms -1 , which occurs at I c± 165° , b ~ 55° , as well as on the opposite side of the sky. 



component as 



(19) 



where Vq is a 3 x 3 matrix. In order for both the divergence and curl of vg(r) to vanish, Vq must be a traceless, 
symmetric matrix. Consequently, it has only five independent elements, two diagonal and three off-diagonal. 

We could allow for the presence of such a quadrupole in velmod by treating these five elements as free 
parameters. However, this is a dangerous procedure, because the modeled quadrupole would then have the 
freedom to fit the quadrupole already present in the IRAS velocity field, which is generated by observed density 
fluctuations. We wish to allow for the external quadrupole, but we do not want it to fit the /^-dependent 
quadrupolar component of the IRAS-predicted velocity field. In other words, we want the external quadrupole 
to be that required for the true value of Pi, which we do not know a priori, rather than the "best fit" value at 
any given /3j. This problem would indeed be very serious if inclusion of the quadrupole made a large difference in 
the derived value of (3j. Fortunately, however, it does not. As we show below, we obtain a maximum likelihood 
value Pi = 0.56 when the quadrupole is not modeled. When we treat all five components of the quadrupole as 
free parameters for each pi, we obtain Pi = 0.47.0 Because the best-fit quadrupole is relatively insensitive to 
Pi, we can estimate the external quadrupole by averaging the fitted values of the five independent components 
obtained for Pi = 0.1, 0.2, . . . , 1.0. In this way, we "project out" the /^-independent part of the quadrupole. In 
our final velmod run, we use this average external quadrupole at each value of Pi. Throughout, we ignore the 
very small effect that this quadrupole might have on the derived IRAS density field. 

In Figure 4, this quadrupole field is plotted on the sky in Galactic coordinates for a distance of 2000 kms -1 . 
The inflow due to the quadrupole, which occurs near the Galactic poles, is of greater amplitude than the 
outflow, which occurs at low Galactic latitude. The quadrupole reaches its maximum amplitude at I ~ 165° , 
b ~ 55° , in the direction of the Ursa Major cluster, as well as on the opposite side of the sky. In § [5], when 
we plot velmod residuals on the sky with and without the quadrupole, the need for the quadrupole field shown 



This value differs from the value of 0.49 quoted in the Abstract because we will not allow the quadrupole to be free parameters at each value of /3j. 
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Fig. 5. The VELMOD likelihood statistic, .Cforw (Eq. |l3| ) , plotted as a function of /3j for the real data. In the left hand plot, an external quadrupolc 
is modeled, as described in the text. In the right hand plot, no external quadrupolc is included in the velocity field. Cubic fits to the likelihood points 
are shown as dotted lines. The minima of the fitted curves, /3 m i n , are the maximum likelihood estimates of f3j in each case. Note the very different 
values of the vertical axes of the two plots; this indicates the large increase in formal likelihood when the quadrupolc is included. 

in Figure 4 will become clear. Indeed, we will show in § |5| that the velmod fit is statistically acceptable only 
when the quadrupole is included. Table |2| tabulates the numerical values of the independent elements of Vq 
that generate this flow. The rms value of this quadrupole over the sky is 3.3%, pleasingly close to the value we 
expect from theoretical considerations (Appendix ^). 

When both the quadrupole and the Local Group random velocity vector are modeled, the radial peculiar 
velocity u(r) that enters into the likelihood analysis (see Eq. ||) is given by 

u(r) = (v mA s(r) + Vqy - w LG ) ■ r . (20) 

We emphasize again that while the three components of the Local Group random velocity wlg are treated as 
free parameters in velmod, the five independent parameters of Vq are not, with the exception of a single run 
we used to obtain and then average their fitted values at each f3j. In the final run, from which we derive the 
estimate of Pi quoted in the Abstract, the quadrupole velocity field shown in Figure 4 was used at each value 
of/?/. 

4.5. Results 

The outcome of applying velmod to the A82+MAT subsample described above is presented in Figure 5. 
The velmod likelihood curves are shown both with and without the external quadrupole included. The formal 
likelihood is vastly improved when the quadrupole is included in the fit: since the likelihood statistic £f or w is 
defined as — 2 In [.P(data|/3)] , the ~20 point reduction in the minimum of £f OTW , minus the five extra degrees 
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Fig. 6 . Breakdown of the VELMOD likelihood statistic among subsamplcs. The left hand panel plots likelihood per point versus 0i for the A82 

and MAT samples individually. The right hand panel plots likelihood per point versus /3j for three different rcdshift intervals containing roughly the 
same number of objects. In each case, the minimum occurs within ~0.1 in f3j of the global minimum in £f rw at f3j — 0.492. 



of freedom when the quadrupole is modeled corresponds to a probability increase of a factor ~ e 



7.5 



2000. 



The improvement in formal likelihood through the addition of the quadrupole is so pronounced that we take 
the maximum likelihood value of (3i from that fit, 0.492 ± 0.068, as our best estimate. However, the maximum 
likelihood estimate of (3j when the quadrupole is neglected, 0.563 ±0.074, differs from our best value at only the 
1-cj level. While the quadrupole is important, it does not qualitatively affect our conclusions about the likely 
value of f3 1. 

We can make several additional tests of the robustness of our results. Figure 6 shows how the likelihoods per 
object break down for fits to different cuts on the sample; see also Table [|. The left hand panel plots £{ orw /N 
versus for the A82 and MAT samples separately, where N = 300 for A82 and N = 538 for MAT. Cubic fits 
to the individual sample likelihoods yield fa = 0.489 ± 0.084 and 0.498 ± 0.107 for A82 and MAT respectively. 
This agreement is remarkable, given that there are only 53 galaxies in common between the two samples. Note 
that the /3-uncertainty is larger for the MAT sample, even though it contains nearly twice as many objects 
as the A82 sample. This is because the MAT objects typically lie at larger distances than do A82 objects, a 
property of the likelihood fit we now illustrate. 

The right hand panel of Figure 6 plots C{ OTW /N versus @i for three subsamples in different redshift ranges. As 
Table Q shows, the agreement in the derived values of (3i is quite good. Changing the specific redshift intervals 
used for this test does not significantly change the results. Note that the ^-resolution decreases as one goes 
to higher redshift, despite the fact that there are nearly equal numbers of objects in each of the three redshift 
bins. This is because the likelihood is sensitive mainly to the fractional distance error in the IRAS prediction. 
Hence, nearby galaxies are more diagnostic of incorrect peculiar velocity predictions, and thus of /3j. 

The fact that £f orw /iV decreases with redshift should not be interpreted as meaning that more distant 
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Fig. 7. Left hand panels: the TF slopes (top) and scatters (bottom), for the A82 and MAT samples, derived from VELMOD as a function of /3j. 

Right hand panels: the amplitude of the Local Group random velocity vector (top) and the velocity noise a v (bottom) derived from VELMOD as a 
function of j3j. All plots correspond to the run in which the quadrupolc was modeled. 

objects are better fit by the velocity model. This decrease instead reflects a property of the velmod 
likelihood implicit in Eq. (|l5|), which shows that the expectation value of Cf OTW /N is ~ 1 + ln(27r) + 



In 0j F + (2.17cr v /[w(l + u'(w))]) , which increases with decreasing cz in general. This effect will be par- 
ticularly pronounced in flat zones (u 1 ~ —1) in the redshift-distance relation, which are found in the Local 
Supercluster, which is why there is a marked difference between £f orw /iV for the A82 and MAT samples (the 
former preferentially populates the Local Supercluster region). 

In Figure 7 we plot for the real data the same quantities plotted for a mock catalog in Figure 3, as well as 
the TF slopes. The slopes are extremely insensitive to /?/. This indicates that the IRAS assigns low and high 
linewidth galaxies nearly same relative distances at all (3i . Significantly, the amplitude of the fitted Local Group 
velocity vector is minimized near the maximum likelihood value of Pi, just as we saw with the mock catalog. 
This indicates once again that the fit attempts to compensate for a poor velocity field at very low and high 
(5i by moving the Local Group. The mock catalogs showed us that the errors on the Cartesian components of 
wlg are of order 50 km s" 1 (Table []]). Thus, the small value of |wlg| obtained from velmod indicates that the 
Yahil et aJ. (1977) transformation to the Local Group barycenter is correct to within ~ 50 kms" 1 , and that the 
Local Group has random velocity ^50 kms -1 relative to the mean peculiar velocity field in its neighborhood. 

The lower right panel of Figure 7 shows that a v increases monotonically with f3j. Its maximum likelihood 
value is 125 kms -1 . This is a remarkably small number, when one considers that it includes the effect not 
only of random velocity noise but also of IRAS prediction error. In particular, if our estimate of the IRAS- 
prediction errors derived from our mock catalog experiments (§ [^), ~84 kms" 1 , are roughly correct, our 
value for a v implies that the true 1-dimensional velocity noise is ,$,100 kms" 1 . This result is consistent with 
past observations that the velocity field is "cold" (cf., Sandage 1986; Brown & Peebles 1987; Burstein 1990; 
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Fig. O. Top panel: The VELMOD likelihood statistic £f orw as a function of /3j for a two runs in which the Local Group velocity vector was forced 

to vanish, and the velocity noise parameter <r v was held fixed at 150 and 250 kms" 1 . Although the formal likelihoods of the fit are worse than that of 
the full fit (compare with Figure ^|). particularly for a v — 250 kms - 1 , the maximum likelihood estimates of 0i arc nearly unchanged. Bottom panel: 
variation in the TF scatters cttf i for the A82 and MAT samples, as a function of (3j for these VELMOD runs. The larger values correspond to the 
<7v — 150 kms - 1 run. Note that the minimum derived TF scatters do not necessarily correspond to maximum likelihood; sec text for further details. 

Groth, Juszkiewicz, & Ostriker 1989; Strauss, Cen, & Ostriker 1993; Strauss, Ostriker, & Cen 1997). Finally, 
the lower left panel demonstrates again what was seen earlier with the mock catalogs (Figure 3), namely, that 
maximizing probability does not correspond to minimizing TF scatter. In large measure, this is because there 
is a tradeoff between the variance due to the velocity noise a v and that due to the TF scatter. As (3j approaches 
1, a v gets steadily larger; <ttf gets correspondingly smaller, despite the fact that the high /3j models are worse 
fits to the TF data. The TF scatters level out or rise only at /?/ ~ 1. 

A final test of robustness involves eliminating the freedom in the velmod fit provided by the parameters wlg 
and a v . One could argue that these parameters are like the quadrupole: they "are what they are," and we 
should not allow them to absorb the fit inaccuracies at the wrong value of /?/. To assess this, we carried out 
two velmod runs in which wlg was assumed to vanish identically In the first run, we fixed the value of a v 
at 150 kms -1 , and in the second at 250 kms -1 . The quadrupole was held fixed at its best fit value; the free 
parameters in this fit were limited to /3j and the three TF parameters for each of the two samples. The results 
of this exercise are shown in Figure 8 and in Table [|. The derived values of f3j differ inconsequentially from our 
best estimate obtained from the full fit. This shows that allowing ourselves the freedom to fit both wlg and cr v 
does not materially affect the derived value of (3j. The formal uncertainties in 0j are much reduced relative to 
the full fit because formerly free parameters have been held fixed. For a v = 150 kms -1 the formal likelihood is 
worse than for the full fit, but only at the ~ 2a level. This reflects the fact that a v = 150 kms -1 and wlg = 
themselves differ by only ~ 1 a from their maximum likelihood values, according to our error estimates from 
§ 3^. However, the formal likelihood for the a v = 250 kms" 1 run is considerably worse (by a factor of ~ 10 -7 ) 
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than for the full fit. This shows that we rule out such a large a v at high significance. 

The bottom panel of Figure 8 shows the fitted values of otf a s a function of j3i for each of the two values 
of a v and for each of the two TF samples. The TF scatters now track the likelihood much better than they 
did in the full fit (bottom panel of Figure 7); with a v fixed, maximizing likelihood is more nearly equivalent to 
minimizing TF scatter. However, they are still not the same thing: likelihood maximization occurs for f3j ~ 0.5, 
whereas TF scatter is minimized at fli ~ 0.6. This is due to the nonlocal nature of the probability distribution 
described by Eq. (|i~l| ) (cf. Figure 1). The likelihood of a given data point depends on the peculiar velocity and 
density fields all along the line of sight interval allowed by the TF and velocity dispersion probability factors, 
not merely on how close the TF-inferred and IRAS-predicted distances are to one another. 

The bottom panel of Figure 8 also shows that the TF scatter one derives from velmod depends on the value 
of a v . The full fit told us that IRAS errors plus true velocity noise amount to ~ 125 kms -1 . The values of cttf 
obtained in the full fit (Table ^) absorbed the remaining variance. Changing a v to 150 kms" 1 reduces the TF 
scatters by about 0.01 mag. With o v fixed at 250 kms" 1 , however, we find 0.39 and 0.40 mag for the A82 and 
MAT TF scatters. While these latter values are certainly underestimates, the large changes demonstrate that it 
is very difficult to estimate cttf to high accuracy because of its covariance, however slight, velocity noise. This 
is one reason that it is inadvisable to use the value of cttf obtained from fitting TF data to peculiar velocity 
models as a measure of the goodness of fit. We return to this issue in § || below. 

4.6. velmod Results using 500 kms 1 Smoothing 

In our mock catalog tests, we found that using 500 kms -1 smoothing in the IRAS reconstruction resulted 
in ~ 20% overestimates of /?/ (§ p.l| ). However, because the mock catalog may not faithfully reproduce the 
dynamics of the real universe, it is useful to see how much /?/ changes for the real data when the 500 kms" 1 - 
smoothed IRAS reconstructions are used. We carried out two such velmod runs, one with and one without the 
quadrupole. (We determined the quadrupole the same way as for the 300 kms -1 smoothed reconstruction, and 
found that the two differ little.) The resulting maximum likelihood estimates of /3j are listed in Table 0. The 
larger smoothing results in an increase in /?/, as expected. However, the 500 kms -1 result, j3i = 0.544 ± 0.071, 
is within 1 a of our favored result obtained at 300 kms" 1 smoothing. If we reduce this value by 20% in accord 
with the bias seen in the mock catalogs, we obtain f3j = 0.45 ± 0.07, also within 1 a of our preferred result. 
Our choice of a 300 kms -1 smoothing scale is thus unlikely to have led us seriously astray, even if the mock 
catalogs are imperfect guides. 

4.7. Consistency of the Mark III and velmod TF Relations 

In constructing the Mark III catalog, Willick et al. (1996, 1997) required that the TF distances for objects 
common to two or more samples agree in the mean. As noted above, velmod yields an independent TF 
calibration for each sample included in the analysis. As a further consistency check, we can ask whether the 
velmod TF calibrations for the A82 and MAT samples are also mutually consistent. 

We compared A82 and MAT TF distances using the velmod TF relations for 75 objects common to the two 
samples. We limited the comparison to objects whose A82 versus MAT TF distance moduli differ by 0.8 mag 
or less. (Not all of these objects were part of the velmod analysis, as some did not meet the criteria outlined 
in § ^^). We found that the velmod calibrations yield an average distance modulus difference (in the sense 
MAT-A82) (Afj) = -0.056 ± 0.046 mag; the Mark III TF calibrations yield (A/z) = 0.018 ± 0.046 mag. The 
corresponding median distance modulus differences are —0.015 mag (velmod) and 0.035 mag (Mark III). Thus, 
as measured by the criterion of generating mutually consistent TF distances among samples, velmod gives the 
correct result. In Table ^ we list the velmod TF parameters and their Mark III counterparts. We see that 
the A82 zero points, slopes, and scatters derived from the two methods are in almost perfect agreement. The 
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Fig. 9.— Same as Figure ^| except that the A82 and MAT TF parameters have been fixed at their Mark III catalog values. The extremely small 
formal error bars result from fixing the TF parameters and arc unrealistic. The maximum likelihood values of 0j differ negligibly from those obtained 
when the TF parameters are free. 

MAT zero points and scatters also agree to well within the errors. The MAT slopes show a somewhat larger 
discrepancy. However, the two slopes are nearly within their mutual error bars; moreover, the MAT sample 
use here is only about half as large as that used by Willick et al. (1996) in deriving the MAT TF slope. In any 
case, this slope difference, even if real, is of no consequence for determination of f3j, as we now show. 

As a final test of VELMOD-Mark III consistency, we ran velmod without allowing the TF parameters to vary, 
instead holding them fixed at their Mark III values. We did so both with and without the quadrupole, while 
holding c v fixed at 150 kins -1 and setting wlg = (note from Figure 8 that these latter velocity parameters 
yield the same f3j as when they are allowed to vary freely). The results of this exercise are shown in Figure 9 
and tabulated in Table |2[ As can be seen, while there is a large formal likelihood decrease relative to the best 
solution, using the Mark III TF relations has a negligible effect on the value of /3j obtained from velmod. In 
particular, use of the Mark III TF relations does not bring our velmod result appreciably closer to the potiras 
result, Pi = 0.86, of Sigad et al. (1997). We discuss this issue further in § |6.1.1 . Note that, in contrast with 



full velmod, neglect of the quadrupole now has no effect on the derived although its inclusion still results 
in a significant likelihood increase. The indicated formal error bars on /3j should not be taken literally here, 
because fixing the TF zero points prevents them from compensating for IRAS zero point errors (cf. 



5. Analysis of the Residuals: Do the Predictions Match the Observations? 

The velmod analysis can tell us which velocity field models — which values of f3j, o~ v , wlg, and quadrupole 
parameters — are "better" than others. However, as with maximum likelihood approaches generally, it cannot 
by itself tell us which, if any, of these models is an acceptable fit to the data. This is because we do not have 
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precise, a priori knowledge of the two sources of variance, the velocity noise a v and the TF scatter <ttf- We 
have instead treated these quantities as free parameters and determined their values by maximizing likelihood. 
As a result, a standard x 2 statistic will be ~ 1 per degree of freedom even if the fit is poor. 

We can, of course, ask whether the values of <7tf and o v obtained from velmod agree with independent 
estimates. It is reassuring that they do. We find cttf — 0.46 mag for both the A82 and MAT samples, 
within the range estimated by Willick et aJ. (1996) by methods independent of peculiar velocity models. This 
agreement is of limited significance, however. TF scatter is very sensitive to non-Gaussian outliers (§ 4.1), and 
thus to precisely which objects have been excluded. Furthermore, the MAT subsample used here is only about 
half as large as the MAT subsample used by Willick et al. (1996) to estimate its scatter. The velmod result 
for the velocity noise, a v ~ 125 kms -1 , is remarkably small, and appears consistent with recent studies based 
on independent methods (e.g., Miller, Davis, & White 1996; Strauss, Ostriker, & Cen 1997). Indeed, because 
~ 90 kms -1 may be attributed to IRAS velocity prediction errors (§ |3.2j ), our value of a v suggests a true 1-D 
velocity noise of ^90 kms" 1 . Still, the small o~ v is not necessarily diagnostic; for demonstrably poor models 
(e.g., fii < 0.2) we find an even smaller value of a v . An alternative approach is thus required for identifying a 
poor fit. 

Consider fitting a straight line y = ax + b by least squares to data (xi,yi) whose errors are unknown. One 
obtains a, b, and also the rms scatter about the fit. Because the scatter is derived from the fit, the x 2 statistic 
is ~ 1 per degree of freedom by construction. However, if the straight line is a bad fit — if, say, the relation 
between y and x is actually quadratic — then the residuals from the fit will exhibit coherence. Coherent residuals 
in excess of what is expected from the observed scatter would signify that a model is a poor fit. In this section, 
we will make such an assessment for the velmod residuals. We will first define a suitable residual and plot it 
on the sky. We will demonstrate coherence and incoherence of the residuals, for "poor" and "good" models 
respectively, by plotting residual autocorrelation functions. Motivated by these considerations, we will define 
and compute a statistic that measures goodness of fit. 



5.1. Sky Maps of velmod Residuals 



velmod does not assign galaxies a unique distance (§ 2.2.1). Thus, there is no unique measure of the amount 



by which their observed and predicted peculiar velocities differ. However, there is a well-defined expected 
apparent magnitude for each object, 



E(m\rj, cz) = / mP(m\n,cz)dm , (21) 

J — oo 

where P(m\r], cz) is given by Eq. ([[l]). Similarly, the rms dispersion about this expected value is 



Am = y E(m 2 \r], cz) — [E(m\r], cz)] . (22) 

Note that Am is not equal to the TF scatter; it also includes the combined effects of velocity noise, peculiar 
velocity gradients, and density changes along the line of sight. At large distances, Am tends toward cttf 
(although dispersion bias can make it smaller; cf. Willick 1994). With the above definitions, one can define a 
normalized magnitude residual for each galaxy 

m- E(m\mcz) 

m = t • (23) 

Am 

The normalized residual has the virtue of having unit variance for all objects. In contrast, the variance of the 
unnormalized magnitude residual m — E(m\n,cz) depends on distance (velocity noise is more important for 
nearby objects), while the variance of a peculiar velocity residual formed from the magnitude residual (Eq. |24| 



25 



Real data residuals, 0<cz^ 1000 km s 1 



360; 




-190- 



-;0 



f]=0.G (No Quad) 



: cjt 



) 150 km s- 1 
^250 km a- 1 



360; 



ft=0,5 (Quad) 




-1?0- 



.jo 



) 150 km s" 
7250 km s~ 



Fig. 10 . VELMOD velocity residuals plotted on the sky in Galactic coordinates, for objects with < cslg £ 1000 kms 1 . The top panel is for 

the j3j — 0.6 run without the quadrupolc. The bottom panel is for the j3j — 0.5 run with the quadrupole modeled. Open circles indicate objects that 
arc inflowing relative to the IRAS prediction; stars indicate outflowing objects. 

below) grows with distance. The normalized magnitude residual 5 m is a measure of the correctness of the 
IRAS velocity model. If, in a given region, 5 m !> in the mean, galaxies in that part of space must be more 
distant than IRAS predicts them to be — i.e., they have negative radial peculiar velocity relative to the IRAS 
prediction. Regions in which 5 m < in the mean have positive radial peculiar velocities relative to IRAS. 

We will use the normalized magnitude residual below in our quantitative analysis of residuals, but let us 
first visualize this residual field on the sky, by converting 8 m into the corresponding radial peculiar velocity 
residual 6u. Were we to do this for each galaxy individually, the ~ 20% distance errors due to TF scatter would 
completely hide systematic departures from the IRAS model. Instead, we will compute smoothed velocity 
residuals. This procedure is most well-behaved if we first smooth S m and then convert the result into 5u. We 
first place each galaxy at the distance d assigned it by the IRAS velocity model; this is a redshift space 
distance so our calculation is unaffected by Malmquist bias. Then for each galaxy i, we compute a smoothed 
residual 8^ ni as the weighted sum of the residuals 5 m of itself and its neighbors j, where the weights are 

Wij = exp (— dij/2Sf J , and dij is the ffiAS-predicted distance between galaxies i and j. We take the smoothing 
length Si to be Si = di/5. The smoothed residual 5 s mi is converted into a smoothed velocity residual according 
to 



5ui = di h -/ao - 2( ^ xAmi) 



(24) 



Wc take d to be the "crossing point distance" w defined in § 2.2.2 In the case of triple- valued zones, we take the central distance. 
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Fig. 11. Same as Figure E3, but for objects with 1000 < cz LG < 2000 km s 1 . 

where Ami is given by Eq. (|22"|). The quantity fi is given by exp (— A?/2) , where Aj = 0.46Amj/ ^Jj2j w ij'i ^ 
guarantees that 5uf, which is log- normally distributed, has expectation value zero if Aim (which is normally 
distributed) does (cf. Willick 1991, § 6.3, for details). 

In Figures 10, 11, and 12 we plot velmod velocity residuals on the sky for the redshift ranges 0-1000 kms -1 , 
1000-2000 kms -1 , and 2000-3000 kms -1 respectively. In each fi gure, the top panel shows residuals from the 
Pi = 0.6 (no quadrupole) fit, and the bottom panel shows residuals from the /?/ = 0.5 (quadrupole modeled) fit, 
the velmod runs closest to the maximum likelihood value of (3i for each case. The plots reveal why the addition 
of the quadrupole results in a large increase of likelihood. In each redshift range, the no-quadrupole fits show 
coherent negative velocity residuals in both the Ursa Major region [l ~ 150° , b ~ 65° ), and at b ~ —60° , 
I ^ 30° and / ^ 330° . In both of these regions, the addition of the quadrupole greatly reduces the amplitude 
of the residuals. In other parts of the sky, smaller but still significant coherent residuals are reduced with the 
addition of the quadrupole. This shows that the pattern of departure from the pure IRAS velocity field is 
well-modeled by a quadrupolar flow of modest amplitude, and therefore has the simple physical interpretation 
we discussed in § |4.4j . 

In the bottom panels, it is difficult to find any well-sampled region within 2000 km s -1 where \Su\ ^ 100 km s -1 . 
This is all the more remarkable because the TF errors themselves are of order 300 kms -1 per galaxy at a 
distance of 1500 kms -1 . Figure 12 does show several high-amplitude residuals. However, at 2500 kms -1 , the 
TF residual for a single object is 500 kms -1 , so when the effective number of galaxies per smoothing length 
is only a few, velocity residuals of several hundred kms -1 are expected from TF scatter only. In well-sampled 
regions, one sees that in general \5u\ ^150 kms -1 , the only exception being a patch of large (^250 kms -1 ) 
positive residuals at / ~ 330° , b ~ —20°. In the b > 0° part of the Great Attractor region at I ~ 300° , the 
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Real data residuals, 2000<cz^3000 km s 





Fig. 12. Same as Figure [uj but for objects with 2000 < cz LG < 3000 km s 1 . 

residuals are < 100 kms -1 even in this highest redshift shell. This is significant, given the oft-heard claims 
that the IRAS model cannot fit the observed flow into the Great Attractor. 

In Figures 13, 14, and 15 we again plot velmod residuals on the sky for the three redshift ranges, now for 
the two values of 0i most strongly disfavored by the likelihood statistic in the range studied, /3j = 0.1 (top 
panels) and (3j = 1.0 (bottom panels). In each plot, the quadrupole of Figure 4 has been included. These 
plots, which should be compared with the bottom panels of Figures 10, 11, and 12, demonstrate why very 
low and high (3j do not fit the TF data well. In each redshift range, these models exhibit large, coherent 
residuals. For Pi = 0.1, we see large negative peculiar velocities relative to IRAS in the Ursa Major region 
at cz < 2000 kms -1 . Indeed, the residual plot for (3j = 0.1 (with quadrupole included) shows many of the 
same features as the no-quadrupole model with (3i = 0.6, because the IRAS field itself contributes some of 
the needed quadrupole. However, the IRAS contribution scales with and is thus inadequate at low fii. At 
(5i = 1.0 many of the systematic residuals associated with the quadrupole are gone, especially in Ursa Major. 
However, other regions show highly significant residuals: at I ~ 150° , b ~ —20° and cz < 1000 kms -1 , for 
example, one sees negative peculiar velocity residuals of amplitude ^200 kms -1 , which is significant at such 
small distances. In the same redshift range, at I = 270-360° , b < 0° there are positive velocity residuals of 
amplitude ^ 150 kms -1 . These regions exhibit much smaller residuals in the Pi = 0.5 model. 

In the higher redshift shells, the poor fit of the /?/ = 1.0 model is evidenced chiefly in the direction of the 
Great Attractor (I ~ 300° , b ~ 20° ). For 1000 < cz < 2000 kms -1 , this model predicts much too large positive 
peculiar velocities, so that the data exhibit inflow relative to the model. In the highest redshift bin, the (3i = 1.0 
model exhibits both positive and negative velocity residuals of high amplitude in the GA direction; residuals 
of both signs are seen in this region for (5i = 0.5 as well, but they are of much smaller amplitude (lower panel 
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Real data residuals, 0<cz^ 1000 km s 1 
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Fig. 13. Same as Figure jicj, except now results for /3j — 0.1 and @i — 1.0 arc shown. In each case, the quadrupolc is the same as it was for the 

best fit model (0i = 0.5). 

of Figure 12). The /3j = 0.1 model, on the other hand, predicts too-small positive peculiar velocities in the 
GA direction at the highest redshifts. Indeed, note that in the 2000 < czlg < 3000 kms" 1 shell, nearly all 
data points exhibit outflow relative to the /?/ = 0.1 IRAS predictions, whereas at lower velocities the residuals 
typically indicate inflow. This global mismatch is more general than the insufficient quadrupole mentioned in 
the previous paragraph, showing that low /3j could not yield a good fit even if we were to give velmod full 
freedom in fitting the quadrupole at all f3j. 

Although sky plots of residuals argue in favor of the f3j = 0.5 plus quadrupole model, the residuals from that 
model are not manifestly negligible. We will address this issue quantitatively below. For now, however, we can 
demonstrate qualitatively that the residuals seen in the /3j = 0.5 plus quadrupole model are not unexpected by 
comparing with the mock catalogs, for which the IRAS velocity predictions are known to be a good fit. Figure 16 
plots velmod velocity residuals with respect to /3j = 1 (the correct value) for a single mock catalog. The same 
three redshift ranges used for the real data are shown. The mock catalog residuals are comparable in amplitude 
and apparent coherence to the real data. Generally speaking, velocity residuals in well-sampled regions are 
£ 100 kms -1 within 1000 kins" 1 , and are ,$200 kms -1 at larger distances. One also sees apparent coherence 
in the mock catalog residual map, as was the case with the real data. The similar amount of apparent coherence 
in the real and mock data indicates that the former is not a result of a poor fit. The apparent coherence in the 
residual sky maps is an artifact of the smoothing used to generate them, as we show in the next section. 
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Real data residuals, 1000<cz^2000 km s 1 
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Fig. 14.— Same as the previous figure, but for objects with 1000 < cslg < 2000 km s 

5.2. Residual Autocorrelation Function 

The sky plots shown above provide visual evidence that the Pi = 0.5 plus quadrupole fit has generally small 
residuals, although they are correlated to some degree. In this section, we quantify these correlations with the 
residual autocorrelation function 



1>{t) = 



1 



E 



3m,i3m,j j 



(25) 



=T±Ar 



where 5 m was defined in Eq. (^), and the sum is over the N p (t) distinct pairs with IRAS predicted separation 
dij within At = 100 kms -1 of a given value r. This definition makes i^{t) insensitive to the values of otf 
and o~ v (because the 8 m ,i are themselves normalized using their maximum likelihood values for each /?/), but 
sensitive to the residual correlations that signal a poor fit. 

In Figure 17, we plot ip(r) versus r for the IRAS plus quadrupole models, with (3j = 0.5, 0.1, and 1.0, as well 
as the Pi = 0.6, no quadrupole model. The error bars are described below. The model that fits best according 
to the likelihood statistic, Pi = 0.5 plus quadrupole, shows no significant residual correlations on any scale. 
The correlation function is everywhere consistent with zero, as we would expect if the IRAS velocity field plus 
the quadrupole is indeed a good fit to the data. Indeed, the absence of residual correlations is the basis for a 
statement made in § 2.2.1, namely, that the individual galaxy probabilities P(m\rj,cz) are independent, and 
thus validates the velmod likelihood statistic £f or w 

The other models shown in Figure 17 all exhibit significant residual correlations. The Pi = 0.6, no quadrupole 
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Real data residuals, 2000<cz^3000 km s 1 
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Fig. 15.- Same as the previous figure, but for objects with 2000 < cslg < 3000 km s 

model has noticeable correlations on small and large scales, as does the Pi = 0.1 plus quadrupole model. Indeed, 
several of the values of ip(r) for (3j = 0.1 are so large that they are off-scale on the plot. The Pi = 1.0 plus 
quadrupole model exhibits strong correlations for r^2000 kms" 1 , although it is well-behaved on large scales. 



5.2.1. Using Residual Correlations to Identify Poor Fits Quantitatively 

In order to compare the observed residual correlations with the results from the mock catalogs, we would like 
to define a single statistic that summarizes the deviation of VK 1 ") from unity. Let us define £(r) = N p (t)i{j(t) 



(cf. Eq. 25). In Appendix C, we discuss the properties of this statistic in greater detail. As we show there, £(r) 
approximates a Gaussian random variable of mean zero and variance N p (r), if indeed the velmod residuals are 
uncorrelated on scale r. (This property was used to compute the error bars on ip{r) above.) To the degree this 
approximation is a good one, the quantity 

(26) 

will be distributed approximately as a x 2 variable with M degrees of freedom, where M is the number of 
separate bins in which £(t) is calculated. In contrast, if the residuals are strongly correlated on any scale r, %| 
will significantly exceed its expected value. 

However, because a single galaxy will appear in many different pairs in the correlation statistic, both within 
and between bins in r, the assumptions made above do not hold rigorously. In Appendix C we explore this issue 
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Mock catalog residuals, $f= 1 




Fig. 16 . VELMOD velocity residuals for a single mock catalog run using @j — 1.0, the true value for the mock catalog. (The particular simulation 

used had a maximum likelihood value of f3j — 0.963.) The three panels show residuals for the three rcdshift ranges used in analyzing the real data. 

further. For now, we appeal to the mock catalogs to assess how closely the quantity %| follows x 2 statistics. 
We computed it for each of the 20 mock catalog runs (§ ^) with j3j = 1. We carried out the calculation to a 
maximum separation of 6400 kms" 1 , in bins of width 200 kms -1 , so that M = 32, and found a mean value 
(x|) = 27.83 ± 1.82, which may be compared with an expected value of 32 for a true x 2 statistic. The rms 
scatter in x\ was 8.15, which is the same as that expected for a true x 2 ■ The difference between the mean and 
expected values is 2.3 a, indicating that x 2 1S n °t exactly a x 2 statistic, for reasons discussed in Appendix C. 
However, because the departure from true x 2 statstics is small, x| is a useful statistic for measuring goodness 
of fit when calibrated against the mock catalogs. 

Before presenting x\ f° r the real data, we consider its variation with f3i for the mock catalogs. In Figure 18 
we plot the average value of %| over the 20 mock catalogs at each value of /3j for which velmod was run. 
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Fig. 17. VELMOD residual autocorrelation functions, ijj{r), plotted for 0j = 0.5 plus quadrupole (best fit model), 0j — 0.6, no quadrupole, 

0i — 0.1 plus quadrupole, and 0j — 1.0 plus quadrupole. In the 0i — 0.1 several points at large r have residuals that arc so large that they fall beyond 
the plot boundaries. 

Although the minimum is at Pi = 1, it is not nearly as sharp as is that of the likelihood as function of Pi (e.g., 
Fig. 2); this statistic does not have the power that the likelihood does for measuring (3j. Indeed, for a single 
realization (the open symbols), the statistic has several local minima. However, it is apparent that a xf value 
much greater than its expected true value of ~ 28 will indicate a poor fit of the model to the data. 

In Figure 19, we plot the statistic %| as a function of Pi for the real data, with and without the quadrupole 
included. The horizontal lines indicate the expected value of %|, and the la and 3a deviations from it. Note 
first that the no- quadrupole model does not provide an acceptable fit for any value of Pi. This is not a conclusion 
we could have reached on the basis of the likelihood analysis alone. When the quadrupole is included, the only 
values of Pi that are unambiguously ruled out are Pi = 0.1, 0.2, and 1.0. The best fit model according to 
velmod, Pi = 0.5 plus quadrupole, also has the smallest value of x|. Given the multiple minima seen for one 
mock realization in Figure 18, this is not necessarily deeply significant. The statistic x| is suitable for identifying 
models that do not fit the data, but does not have the power of the likelihood statistic for discriminating among 
those models that do fit. 

In summary: the velmod likelihood maximization procedure is the proper one for determining which value 
of Pi is better than others, but it cannot identify poor fits to our model. The residual correlation statistic %| 
can identify unacceptable fits, but does not have the power to determine which of the acceptable fits is best. 
We have found that the IRAS velocity field with Pi = 0.5, plus the external quadrupole, is both the best fit of 
those considered, and is also an acceptable fit. Values of Pi > 0.9 and Pi < 0.3 are strongly ruled out. 
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Fig. 18 . The residual correlation statistic x$ > defined by Eq. Pq , plotted as a function of (5j for the mock catalogs. The solid symbols show an 

average over 20 mock catalogs; the open symbols show the values obtained for a single mock catalog. 

6. Discussion 
6.1. What is the value of /?/? 

velmod recovers the correct answer, f3j = 1, to < 10% accuracy when applied to the mock catalogs. At (3i = 1, 
the velocity field in the mock Virgo region is significantly triple- valued. Thus velmod, despite being close in 
spirit to Method II, properly treats triple- valuedness. If the strong triple- valuedness one sees at /3j = 1 were 
present in the real universe, velmod would not assign it unduly small likelihood. Nonetheless, when velmod is 
applied to the real universe, it returns a value of fli = 0.492 ± 0.068 (quadrupole modeled). This value is quite 
insensitive to two other quantities treated as free parameters in the velocity field model, the Local Group random 
velocity wlg and the small-scale velocity dispersion a v (§ |4.5| ). Tests with the mock catalogs demonstrated 
that we obtain an unbiased /?/ using a 300 kms _1 -smoothed IRAS reconstruction (§ |3.l[) . However, we found 
that changing to a 500 km s~ ^smoothed reconstruction makes relatively little difference in (5j (§ |4.6| ). Finally, 
neglecting the quadrupole causes Pi to change by only ~ 1 a. Our conclusion that /?/ ~ 0.5 ±0.07 is thus robust 
against systematic effects internal to our method. 

The velmod result is consistent with the relatively low estimates of f3i obtained from the Method II analyses 
of Hudson (1994), Roth (1994), Shaya et ai. (1995)0, DNW, and Schlegel (1996), as well as those derived from 
comparisons of the IRAS density field with the motion of the Local Group (Strauss et ai. 1992b) and from 
some analyses of the redshift-space anisotropy of the IRAS density field (e.g., Hamilton 1993, 1996; Fisher et 



The Hudson and Shaya et ai. papers actually derive /3 op t , which must be multiplied by ~ 1.3 to obtain an equivalent /3/ ; cf. footnote 2. 
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Fig. 19 . The residual autocorrelation statistic defined by Eq. E6l plotted as a function of /3j for the real data, with and without the quadrupole 

modeled. The heavy solid line shows the expected value of the statistic, which was determined by averaging the derived value for 20 mock catalogs. The 
two dashed lines show 1- and 3-cr deviations from this value. Note that when the quadrupole is not modeled, highly significant residual correlations are 
detected for all values of /?/. (The no-quadrupole points for pi — 0.1 and 0.2 are not shown because their values are too large.) 

al. 1994; Cole, Fisher, Sz Weinberg 1995; Fisher & Nusser 1996). However, it is apparently inconsistent with 
estimates of j3i near unity, as have been found by the potiras analysis (Sigad et al. 1997), measurements of 
the potent fluctuation amplitude (Kolatt & Dekel 1997, Zaroubi et al. 1997), and redshift-space distortions 
of spherical harmonic expansions of the density field (Fisher, Scharf, & Lahav 1994c; Fisher 1994; Heavens & 
Taylor 1995). 

6.1.1. Why do velmod and potiras yield different values of (3i ? 

We do not yet have a satisfactory explanation of why velmod and standard Method II analyses characteris- 
tically yield smaller values of f3j than the Method I potiras approach. One possibility is that the differences 
stem from the Method I/Method II distinction. However, velmod corrects the principal drawback of Method 
II, the inability to deal with multivalued or flat zones in the redshift-distance relation. Thus, if the Method 
I/Method II distinction is at the root of the discrepancy, the reason must be more subtle than the drawbacks of 
standard Method II. Sigad et al. (1997) test for biases in potiras using the same mock catalogs as this paper; 
they too find their determination of (3i to be essentially unbiased. The problem could lie with the Malmquist 
bias corrections that are so crucial to Method I (cf. the discussion in Willick et al. 1997). If these corrections 
are underestimated for any reason — e.g., the TF scatter is larger than estimated, or the density fluctuations are 
larger than modeled — a Method I approach will produce too-strong velocity gradients and thus overestimate 
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(3j. However, the TF scatters used by Sigad et ai. (1997) are consistent with those obtained in this paper, and 
the large potent smoothing limits the effect of Malmquist bias in any case. It is thus unlikely that improper 
Malmquist bias corrections strongly affect the value of /3j obtained from potiras. 

An important difference between velmod and potent is the Gaussian smoothing scales employed, 300 and 
1200 kms -1 respectively. These very different smoothings could result in different values of f3j if the effective 
bias parameters on these scales are different. In order to reconcile velmod and potiras, we would need the 
effective bias parameter to decrease by a factor of 1.7 between scales of 300 and 1200 kms -1 . Such a scale- 
dependent biasing has been suggested by the galaxy formation models of Kauffman et al. (1996), but Weinberg 
(1995) and Jenkins et al. (1996) do not find these trends. A recent analysis by Nusser & Dekel (1997) using a 
spherical harmonic expansion of the velocity field finds /?/ = 1.0 for 1200 kms -1 smoothing, but only 0.6 for 
600 kms -1 smoothing, approaching the value we have found in this paper. Such a change of /?/ with smoothing 
scale could signal scale-dependent biasing. 

Still another difference is the volume considered. We have restricted this analysis to cz < 3000 kms -1 (§ ^), 
whereas the analysis of Sigad et al. (1997) extends to 6000 kms -1 ; only ~ 1/3 of the points used fall within 
3000 kms -1 . If, for whatever reason, bj differed locally from its global value, the velmod result could be biased 
low. In a future paper we will extend the velmod analysis to larger distances; however, our preliminary results 
do not show an increase in /?/ when we do so. In addition to probing a larger volume, the Sigad et al. analysis 
uses the full Mark III sample, ellipticals included; the possibility of systematic differences between the TF 
subset we have used in this paper, and the full sample, is difficult to rule out. Finally, it is conceivable that 
the requirement of pre-calibrating TF relations (potent), as opposed to calibrating them simultaneously with 
fitting the velocity field (velmod and Method II generally) accounts for part of the discrepancy. However, fixing 
the velmod TF parameters at their Mark III values has essentially no effect on the derived value of (3j (§ |4.7|) . 
This strongly argues against the notion that a major difference between velmod and potiras is the TF relations 
themselves. 



6.1.2. The effect of cosmic scatter 

The sphere out to 3000 kms -1 is small; the rms value of density fluctuations within spheres of this radius 
is 20% for COBE-normalized CDM. However, this does not propagate to a cosmic scatter error on our derived 
Pi, for two reasons. First, the IRAS velocity field is determined within a sphere of radius 12,800 kms -1 , within 
which the rms fluctuations are only a few percent. Thus the peculiar velocity field is subject to very little 
cosmic scatter. Second, this scatter primarily manifests itself as a monopole term (cf. the discussion in § 4.4), 
and therefore is fully absorbed into the zero-points of the TF relations (§ |3.3| ), having no effect on the derived 
value of (3j. 



6.2. Do the IRAS and TF Velocity Fields Agree? 

An important conclusion of this paper is that the agreement between the predicted and observed peculiar 
velocity fields is satisfactory (§ ||), as it must be if the resulting estimate of /?/ is to be believed. This agree- 
ment is consistent with the hypothesis that gravitational instability theory correctly describes the relationship 
between the peculiar velocity and mass density fields. It also suggests that the linear biasing model, Eq. (|2|), 
is a reasonable description of the relative distribution of IRAS galaxies and all gravitating matter Gaussian 
smoothed at 300 kms" 1 . 
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6.2.1. Comparison with Davis, Nusser, & Willick (1996) 

DNW reached a different conclusion. Comparing the IRAS and TF velocity fields with a Method II approach, 
(itf; cf. § |2.1| ) DNW found that the fields do not agree at a statistically acceptable level. In particular, a x 2 
statistic resulting from a mode-by-mode comparison of the IRAS and itf velocity fields was found to be 100 for 
55 degrees of freedom ,p| DNW argued that the excessive value of their \ 2 statistic resulted primarily from a 
dipole in the TF velocity field that grows with scale, a feature not seen in the IRAS predictions. They cautioned 
that, as a result, their maximum likelihood value of (3j ~ 0.5 was not necessarily meaningful. 

Why do we find agreement between the TF and IRAS data, while the itf analysis of DNW did not? We 
cannot answer this question with assurance, but we can suggest two likely causes of the discrepancy. First, the 
itf analysis requires that the raw magnitude and velocity width data of the different samples be placed on a 
single, uniform system. This was achieved by applying linear transformations to the magnitudes and widths 
of each sample (Willick et al. 1997). Such a procedure in effect links together the TF zero points of samples 
that probe different volumes. Any systematic error in matching the data sets will manifest itself in spurious 
large-scale motions; in particular, the scale-dependent, dipolar flow found by DNW (see, for example, their 
Figures 12 and 13) is fully degenerate with a zero-point error in the relative TF calibrations of Southern and 
Northern sky samples. Second, DNW extended their itf analysis to 6000 kms -1 , whereas we have restricted 
our analysis to czlg < 3000 kms -1 . In so doing, they (like potiras) incorporated several Mark III TF samples 
(W91CL, HMCL, W91PP, CF) not included in the velmod analysis. It is possible that IRAS and Mark III agree 
locally, but progressively disagree at larger distances. Alternatively, the possible zero point errors mentioned 
above could affect mainly those Mark III samples used by DNW but not included here, given the agreement 
we found between the MAT and A82 distances with the velmod calibrations (§ |4.7| ). 

Since we believe that the DNW discrepancy between the IRAS and TF velocity fields may well be a result 
of systematic errors incurred in matching data sets, an effect to which velmod is insensitive, we are inclined to 
give more weight to our present conclusion that the IRAS— TF agreement is satisfactory. However, if in fact the 
matching of data sets by DNW is validated by ongoing observations aimed at providing reliable North-South 
homogenization (cf. Strauss 1996b), it will be difficult to escape their conclusion that the predicted and observed 
velocity fields do not agree on large scales. In that case, it will be necessary to reexamine the conclusions of 
this paper with regard to the value of /?/. 

6. 2. 2. The Role of the Quadrupole 

Our conclusion that the predicted and observed velocity fields agree also depends on the validity of our 
adopted external quadrupole. Figure 19 shows that only with the quadrupole does our goodness of fit statistic 
x| take on acceptable values. We argue in Appendix || that the 3.3% residual quadrupole we see is mostly 
due to the systematic difference between the true and Wiener-filtered IRAS density field on large scales. The 
residual quadrupole in the mock catalogs is appreciably smaller, < 1%, but this can be understood in terms of 
the different amount of power on intermediate scales {2Tr/k ~ 100 h~ 1 Mpc) in the mock catalog and the real 
universe. Thus, the presence of the quadrupole residual is not evidence for a breakdown of our assumptions of 
gravitational instability theory and linear biasing. 

6.3. What is the value of ft? 

Measuring /3 = O 0,6 /6 is an important objective of velocity analysis. Of course, the more important objective 
is determination of f2 itself. There are, broadly speaking, two ways to proceed. 



Note that unlike this paper, DNW assumed a TF scatter a priori, which allows them to define a goodness of fit directly from their x 2 j cr "- the discussion in § 2.2.1 
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6.3.1. Nonlinear Analysis 

One may attempt to break the degeneracy between Q and biasing by extending gravitational instability 
theory to the nonlinear dynamical regime. In an earlier phase of the velmod project, we attempted to do this; 
very preliminary results of this effort were described in SW, § 8.1.2. In brief, the IRAS reconstruction was done 
as described in Appendix A, but a nonlinear generalization of Eq. ([!]), 

vM _ /(«) f ,3, (l + a 2 vW g (r');b} + av (r' - r) 

V(rj " 4 7r J l + ad[6 g (r');b] | r '- r | 3 ' [ ' 

was used to derive peculiar velocities from the redshift survey density field 6 g . In Eq. a = 0.28 and v = (5 2 ) 
is the mean square value of 5 (Ganon et aJ. 1995; cf., Nusser et aJ. 1991). Note that the mass fluctuation 5 
is written as a generic function of 5 g and b, rather than simply as 5 g /b. This is because once we generalize to 
nonlinear dynamics, we must allow for the possibility of nonlinear biasing as well. There are many ways one 
might imagine doing this (SW, § 2.5; Fry & Gaztahaga 1993). Generically, however, all these complications can 
be expanded to second order to yield a correction to Eq. (||): 
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where 7 parameterizes the combined effects of nonlinear dynamics and nonlinear biasing. 

We carried out a suite of velmod runs using predicted peculiar velocities based on Eq. (^) for a range of 
values of 7, both positive and negative. Our hope was that the velmod likelihood statistic would be significantly 
lower for some value of 7 than for the pure linear case. However, to our surprise, we found that the linear 
dynamics/linear biasing reconstruction (7 = 0) gives the best likelihood of all. We are not certain as to why 
this is. Nonlinear dynamics must enter to some degree, because we know for a fact that 5 g is not everywhere 
■C 1, and indeed can be quite large with our small smoothing. (We of course do not know whether nonlinear 
biasing is important.) Nevertheless, the small scatter between the true and IRAS predicted peculiar velocity 
fields for the mock catalogs (§ 3.2) confirms that the linear IRAS velocity field, smoothed on a 300 kms -1 scale, 
is a good match to actual peculiar velocities that arise from gravitational instability, at least in an iV-body 
simulation. 

A possible explanation of this seeming contradiction is as follows. Our method for predicting peculiar veloc- 
ities (Appendix A) entails assigning a smooth, continuous density field from discrete redshift survey data — a 
procedure which takes into account the probability distribution of distance given redshift (Eq. [A^), smooths the 
data with a 300 kms -1 Gaussian, and applies a Wiener filter — and thus reduces small-scale density enhance- 
ments. In doing so, this procedure mimics qualitatively the effects of nonlinear corrections to the velocity-density 
relation. The good match between the IRAS predictions and the actual peculiar velocities suggests that this 
mimicry is in fact fortuitously good, to the degree that formal nonlinear corrections are unnecessary. 



6. 3. 2. Constraining VL from Independent Estimates of 67 

The second way to estimate fi, given our measurement of (3j, is to constrain bj using independent information. 
If biasing is independent of scale (cf. the discussion in § |6.1.1| ), then bj is the ratio of the rms fluctuations of 
IRAS galaxies on an 8fo -1 Mpc scale, a$(IRAS), to the corresponding mass density fluctuations, erg- Fisher et 
al. (1994a) found that a$(IRAS) = 0.69 ± 0.04 in real space. It follows that /5/ can be viewed as a prediction 
of as for a given value of £1 : 

(0.69 ± 0.04) 

a 8 = ^06 ■ ( 29 ) 
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Fig. 20.— Predictions of ctq as a function of O for open (left panel) and flat (right panel) universes. The lines sloping up and to the right show the 
COBE-normalizcd CDM values of <r& as a function of H, for four values of the Hubble constant, 55, 65, 75, and 85 kms -1 Mpc~ 1 , under the assumption 
of a scale-invariant primordial power spectrum. The line sloping down and to the right shows the result from this paper, as — 0.69/3//f2 ' 6 . The shaded 
region represents the l-cr uncertainty in our value of /3j. 

An entirely independent (though highly model- dependent) way to predict as as a function of f2 is to use COBE- 
normalized power spectra for a range of cosmological parameters. Liddle et al. (1995, 1996) have presented 
fitting functions that provide the normalization of CDM power spectra, in open and flat cosmologies, as a 
function of fi, Q\, the Hubble parameter h = Hq/(W0 kms -1 Mpc -1 ), and the primordial power spectrum 
index n, based on the four- year COBE observations (Bennett et al. 1996; Gorski et al. 1996). Eke, Cole, & 
Frenk (1996) used these fitting functions to obtain as by direct integration of the Liddle et al. power spectra, 
and have kindly provided us with their code for doing this calculation. We may thus constrain as by comparing 
the velmod and COBE/CDM predictions of its value, and requiring that they agree to within the errors. This 
will be the case only for a limited range of f2 (the "concordance range"). We emphasize, however, that the 
discusion to follow depends on two uncertain assumptions: first, that the CMB fluctuations measured by COBE 
can be reliably extrapolated down to 8/i _1 Mpc scales; and second, that the bias parameter is scale-independent 
from 3/j- 1 to 8/i _1 Mpc. 

In Figure 20, we compare the two constraints on as for a scale-invariant (n = 1) power spectrum. The 
left hand panel shows results for an open (i.e., A = 0) universe, and the right hand panel for a spatially flat 
(fi+QA = 1) universe. The COBE/CDM predictions (solid lines labeled with the values of the Hubble constant) 
and the constraint from Eq. (^) (shaded region) scale very differently with $7, so that the two together give 
strong constraints on as and thus 0. The shaded region represents the combined velmod error on (5i and the 
error in as(IRAS) from Fisher et al. (1994a). We do not show corresponding error regions for the COBE/CDM 
predictions which result from uncertainty in the COBE normalization, because the error in the predicted as is 
in fact dominated by the allowed range of Hq, which we take to be 55 < Hq < 85 kms -1 Mpc -1 based on a 
number of recent measurements (Sandage et al. 1996; Freedman 1996; Riess, Press, & Kirshner 1996; Mould et 
al. 1996; Tonry et al. 1997; Kundic et al. 1996). 

Figure 20 gives the following constraints for n = 1. For an open model, the concordance range is f2 = 0.28- 
0.46 with the low (high) value corresponding to the highest (lowest) value of Hq considered. For the flat model, 
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it is = 0.16-0.34. Expressed in terms of the IRAS bias parameter, these ranges correspond to bj = 0.92-1.38 
(open) and bi = 0.68-1.11 (flat). We also considered n ^ 1 flat models. For example, with n = 0.9 the 
concordance ranges are f2 = 0.19-0.40, and 0.21-0.45, depending respectively on whether tensor fluctuations 
are not, or are, included in the COBE normalization (Liddle et al. 1996b). The corresponding bias parameters 
are bj = 0.74-1.21 and bj = 0.80-1.29. 

Two salient points follow from this comparison. First, if Hq > 60 kms -1 Mpc -1 , the concordance range 
for the flat, n = 1 models requires f2,$,0.30, implying U\^0.70. However, studies of gravitational lensing 
have placed an upper limit of S7a < 0.65 at 95% confidence (Maoz & Rix 1993; Kochanek 1996), while a 
recent analysis of intermediate-redshift Type la Supernovae (Perlmutter et al. 1996) indicates < 0.50 
at 95% confidence (both of these constraints apply when a flat universe is assumed). This contradiction 
consitutes evidence against a flat universe with a scale-invariant primordial power spectrum index and Hq > 
60 kms -1 Mpc -1 . If n < 1.0, one can more easily accommodate flat universes with Sl\ < 0.65, provided 
the Hubble constant is ^70 kms" 1 Mpc -1 . The second point is that the combined velmod and COBE/CDM 
predictions of as are extremely difficult to reconcile with an Einstein-de Sitter universe for most reasonable values 
of the remaining cosmological parameters. If one assumes n > 0.9, a Hubble constant ^30 kms -1 Mpc -1 , far 
below current observational limits, would be required for the concordance range to include fi = 1. Alternatively, 
if H = 50 kms" 1 Mpc" 1 , one would require a primordial power spectrum index n = 0.7 and tensor fluctuation 
contributions to the CMB anisotropies. Such a power-spectrum index is at the lowest end of the range currently 
considered plausible in inflationary universe scenarios (e.g., Steinhardt 1996). 

6.4. Summary 

We have described a new maximum likelihood method, velmod, for comparing Tully-Fisher data with pre- 
dicted peculiar velocity fields from redshift surveys. We implemented the method for a czlg < 3000 kms -1 TF 
subsample from the Mark III catalog (Willick et al. 1997), and velocity fields predicted from the 1.2 Jy IRAS 
redshift survey (Fisher et al. 1995). The velocity field prediction is dependent on the value of (3j = Q°' 6 /bj, 
where bi is the bias parameter for IRAS galaxies at 300 kms" 1 Gaussian smoothing. We maximized likelihood 
with respect to the parameters of the TF relation, and several other velocity parameters. 

We applied our method to 20 mock Mark III and IRAS catalogs constructed to mimic the properties of the 
real data. The mock catalogs were drawn from an = 1 iV-body simulation and were constructed to ensure 
bi = 1. Thus, the mock catalogs satisfy (3j = 1. Our velmod runs with the twenty mock catalogs returned a 
mean value of (3i = 0.984 ± 0.018, consistent with the statement that velmod yields an unbiased value of /?/. 
In addition, our mock catalog tests enabled us to assign reliable 1-cr errors to our estimates of (3i, and showed 
that our other derived parameters, including those of the TF relation and the small-scale velocity noise, are 
also unbiased. Because the mock catalogs came from an Q = 1 universe, triple valued zones in the mock Virgo 
region were strong, but were properly handled by the velmod analysis. 

When velmod was applied to the real Mark III data, a considerably smaller value of /3j was derived. If we 
assume that the IRAS-predicted velocity field fully describes the actual one, we obtain (3j = 0.563 ± 0.074. 
However, the residuals from this fit were large and coherent; fitting them by a quadrupolar flow gave a maximum 
likelihood value of [3[ = 0.492 ± 0.068. The quadrupole points toward the Ursa Major cluster, and has an rms 
amplitude of 3.3% of the Hubble flow. In Appendix B, we show analytically that a quadrupole of this amplitude 
is expected given the way that we smooth the density field; its presence is not a sign that the IRAS galaxies 
do not trace the mass responsible for the local flow field. An analysis of the fit residuals demonstrated that the 
IRAS-predicted peculiar velocity field, with the external quadrupole, provides a statistically acceptable fit to 
the TF data within 3000 kms" 1 . The data are thus consistent with the hypothesis that the peculiar velocities 
are due to the gravitational effects of a mass distribution that is proportional to the IRAS galaxy distribution. 
We also find that the data are consistent with a very quiet flow field; the one-dimensional rms noise in the 
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velocity field relative to the IRAS model is 125 ± 20 kms -1 . 

The value of /3j obtained here may also be thought of as a measurement of the rms mass density fluctuations 
ag as a function of f2. Similarly, COBE- normalized CDM power spectra predict a value of as as a function of O 
and other cosmological parameters. If we require that the velmod and COBE- normalized calculations agree, we 
can constrain the value of O. For scale invariant, A = universes, we derive the constraints 0.28 0^:0.46 for 
for 85 ^ Hq 55 kms -1 Mpc -1 . For scale-invariant, flat universes we find 0.16 ^,fl^ 0.34 for the same range of 
Hq. The constraints on Q, shift to higher values (§ |6.3.2 ) if the primordial power spectra are "tilted," n < 1, and 



if tensor fluctuations are present. However, both extreme tilt (n < 0.7) and a Hubble constant at the lowest 
end of the observationally allowed range (Hq < 50 kms -1 Mpc -1 ) would be required to reconcile these results 
with an Einstein-de Sitter universe. 

The conclusions of the previous paragraph all rest, of course, on the validity of our measurement of /?/. 
Tests with mock catalogs show that, subject to our basic assumptions, this measurement is reliable to within 
the quoted errors. We have identified two ways these assumptions can break down. First, the effective bias 
factor bi could depend on scale. In that case, our measurement of /?/, which reflects a 300 kms -1 Gaussian 
smoothing scale, might not be the same as a measurement obtained at larger smoothing; it would not then be 
valid to equate the estimate of as obtained from Eq. |2^ with the COBE/CDM prediction. Second, although we 
have found agreement between the predicted and observed peculiar velocities within 3000 kms -1 , DNW found 
disagreement on larger scales. If the DNW result is validated by future observations (Strauss 1996b) aimed at 
improved TF calibration across the sky, our present claim of TF-IRAS agreement will be undermined. 

There are several areas for further work. One, alluded to in several places in this paper, is to extend our 
analysis to larger redshift, using both the forward and inverse forms of the TF relation. This can be done both 
the Mark III data, and with the extensive new TF (Mathewson et al. 1994; Giovanelli et al. 1997) and D n -a 
(Saglia et al. 1996) samples that are being compiled. We should also consider extending this work to other 
distance indicators; surface brightness fluctuation galaxies (Tonry et al. 1997), with their accurate sampling of 
the nearby velocity field, are natural candidates for the velmod analysis. On the modeling side, this work has 
left us with several conundrums, the most puzzling of which is why the linear IRAS model does so well with 
a smoothing scale of 300 kms -1 . More work is needed with iV-body simulations to understand this. Finally, 
we will not have a coherent picture of the relationship between the velocity and density fields until we can 
understand the different values of f3j obtained by velmod and potiras. 

We thank Marc Davis, Carlos Frenk, and Amos Yahil for extensive discussions of various aspects of this 
project, as well as the support of the entire Mark III team: David Burstein, Stephane Courteau, and Sandra 
Faber. JAW and MAS are grateful for the hospitality of the Hebrew University in Jerusalem, Lick Observatory 
at the University of California, Santa Cruz, and the Astronomy Department of the University of Tokyo for visits 
while we worked on this paper. MAS gratefully acknowledges the support of an Alfred P. Sloan Foundation 
Fellowship. This work was supported in part by the US National Science Foundation grant PHY-91-06678, 
the US-Israel Binational Science Foundation grants 92-00355 and 95-00330, and the Israel Science Foundation 
grant 950/95. 



A. The IRAS Velocity-Density Reconstruction 

The redshifts of galaxies in the IRAS sample are affected by the same peculiar velocities that one is attempting 
to measure in the Mark III dataset. If we measure redshifts cz in the rest frame of the Local Group, then: 

cz = r + r • [v(r) - v(0)] , (Al) 

where v(0) is the peculiar velocity of the Local Group, and v(r) is the peculiar velocity at position r. Indeed, 
because the galaxy density field shows coherence, the galaxy density field measured in redshift space 5 g (s) differs 
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systematically from that in real space, S g (r), as was first described in detail by Kaiser (1987; cf., SW; Strauss 
1996a for reviews). Linear perturbation theory assuming gravitational instability enable us to correct for the 
effects of these velocities. We use here the iteration technique described by Yahil et al. (1991) and Strauss et 
al. (1992c), as updated by Sigad et al. (1997). The density and velocity field are calculated within a sphere 
of radius 12,800 kms -1 ; the density fluctuation field is assumed to be zero beyond this radius. Here we very 
briefly reiterate the improvements described in the Sigad et al. paper, and emphasize certain differences from 
the approach there. 

In regions in which the IRAS velocity field model predicts a non-monotonic relation between redshift and 
distance along a given line of sight, it becomes ambiguous how to assign a distance to a galaxy given its redshift 
(Figure 1). Our approach is similar to that used throughout this paper: we use our assumed density and 
velocity field to calculate a probability distribution of a galaxy along a given line of sight. 

Along a given line of sight, we ask for the joint probability distribution of observing a galaxy along a given 
line of sight, with redshift cz, flux density / and (unknown) distance r: 

P(cz, f, r) = P(cz\r) x P(f\r) x P(r) ; (A2) 

compare with Eq. (|5|). The first term is given by our velocity field model along the line of sight, and is thus 
given by Eq. (||). For the iteration code, we set cr v = 150 kms -1 , independent of position, similar to the best 
fit value we find when we fit for a v from the velocity field data. 
The second term is given by the luminosity function of galaxies: 

P(f\r) = ^(L = A7rr 2 uf)^^rH(L) , (A3) 

df 

where the derivative is needed because the probability density is defined in terms of /, not LF| Finally, the 
third term in Eq. (A2) is given by the galaxy density distribution along the line of sight, Eq. (|8j). 



As described in Sigad et al. (1997), the calculations of the velocity and density fields are done on a Cartesian 
grid. Our approach therefore is to assign each galaxy to the grid via cloud-in-cell (weighting by the selection 
function, of course), where (unlike Sigad et al. 1997) we distribute each galaxy along the line of sight according 



to the distribution function of expected distance, Eq. A2. In order to calculate the selection function for an 



object, we of course need to have a definite position for it; for this purpose, we assign it the expectation value 
of its distance, following Sigad et al. (1997): 

JrP(czJ,r)dr 

W jP(cz,f,r)dr { ] 

Sigad et al. (1997) discuss the use of various filtering techniques to suppress the shot noise in the derived 
density and velocity fields. While they argue for the use of a power-preserving filter for the comparison of the 
IRAS and potent density fields, we have found through extensive experimentation with mock catalogs that 
for the velmod analysis, a Wiener filter gives the best comparison between the density field and the peculiar 
velocity data. 

Finally, we found that when the iteration technique was run to values of fiitl, the density field became 
unstable in the regions around triple-valued zones, oscillating between iterations. We were able to suppress 
these by averaging the derived density field at each iteration with that of the iteration preceding it. This has 
no strong effect on the derived density field for (3 < 1. 



J Eq. (144) of Strauss & Willick (1995) mistakenly left off this last term. 
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B. The Residual Quadrupole 

In this Appendix, we calculate the expected amplitude of the velocity quadrupole generated by density 
fluctuations both external to the IRAS sample (i.e., outside of R = 12, 800 km s -1 ), and internal to it, due to 
the difference between the true density field and the noisy, smoothed estimation of the density field we have 
from the IRAS redshift survey. The IRAS excluded zone is another potential source of quadrupole error, but 
it is filled in by interpolation from regions above and below the excluded zone (Yahil et al. 1991), a procedure 
which agrees well with a multipole interpolation procedure based on spherical harmonics, at least for the 10° 
wide IRAS zone of avoidance (Lahav et al. 1994). 

B.l. The Quadrupole Induced by Fluctuations Beyond the IRAS Volume 

We express peculiar velocity in terms of a potential function <3?(r), such that the radial component of the 
velocity field is given by u(r) = —d^/dr. We will isolate the quadrupole component of this potential, and 
calculate its angle- averaged rms contribution. 

The contribution to from material at distances > R is given by 

m) f ^ 8{v') 

\r'\>R 



$(r) = -^/ d 6 r'-^-. (Bl) 
4vr 7ir'l>R r-r' 



Here, 5 is the mass, not the galaxy, density fluctuation. We now expand the denominator in the integrand in 
terms of spherical harmonics (e.g., Jackson 1976, Eq. 3.70) and isolate the quadrupole term to obtain 

^2 2 j i 



*«(r) = E IW«) f"* I ^ S(r')Yi m (u') . (B2) 

D m=-2 JR r J 



where uj is solid angle. Taking the radial component of the quadrupole velocity uq = —d&Q/dr, squaring, and 
averaging over solid angle gives, after several steps of algebra: 



7/ 2 



= (^f^) 2 E C ^ C ^ J dcoY 2m (u)Y 2 * m ,(") = (^f^) 2 ^El^| 2 ' (B3) 

mm' ' m 

where the last step follows from the orthonormality of the Yi^s, and for convenience we have defined the five 
complex coefficients 



r°° dr' r 

C 2m [R; 5]= — / dJ 5(r') F 2 * m (o/) . (B4) 

J Ft ^ J 



i-°°dr' 
Ir r' 

The expectation value of | C 2rri | 2 is independent of m, so when we take the expectation value of Eq. (|B3|) , we 
can replace the sum with 5 times (C 20 ): 

2 



a . /2/(tt)rV 5 



«>,rm»> = ^-^^J X ■ (B5) 

Using the definition of Y 2 q in terms of the second Legendre polynomial P 2 in Eq. (B4) and Eq. (B5) gives 

, poo A r poo A r pi pi 

< 2 Q,rms(r)) = (m)r) 2 11 ±1 d^d^PMPM^-nl). (B6) 

J R JR ^2 J — 1 J — 1 
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Expressing the correlation function £ as the Fourier Transform of the power spectrum P(k) (e.g., SW, Eq. 46) 
allows the integrals over ri and r 2 to separate. This yields 

<<4,rn»> = U ^ly J d 3 kP(k) W 2 {kR) , (B7) 

where the kernel is given by 

W(kR) =r±r d ^ pM = _ 2 r dr h<M = , (B8) 

Jr r J-i Jr r kR 



and j n is the n-th order spherical Bessel function. Comparison of Eqs. (B7) and flB8| ) with Eqs. (37) and (38) 
of SW allows us to recast our result as 

r- l U rms (r)) 1/2 = 2 -^a R (B9) 



for the expected rms quadrupole velocity on a sphere due to mass density fluctuations at distances > R, 
expressed as a fraction of Hubble flow. Here a\ is the variance in the mass overdensity within spheres of 
radius R. As mentioned in the text, this gives a fractional quadrupole of the order of 1-2% for a variety of 
COBE-normalized power spectra. 



B.2. The Effects of Wiener Filtering and Shot Noise 

The Wiener filter operates on the Fourier Transform of the IRAS density field. The final density field differs 
from the true density fiels for two reasons: the discreteness of the galaxy distribution gives rise to shot noise, 
and the Wiener filter, while suppressing shot noise, also suppresses the density field itself. We calculate the 
contribution to the quadrupole from both effects. 

Let 5r(k) represent the true Fourier component of the underlying (noiseless) density field at wavevector k; 
the quantity with which we calculate the velocity field is the Wiener-filtered noisy image, whose Fourier modes 
are given by: 

J(k) = h(k,r)[5 T (k) + e(k)], (BIO) 
where the Wiener filter itself is (e.g., Zaroubi et al. 1995): 

ft(t - r > = p W+ ^V))- - (Bn » 

and P(k) is set a priori; we used a functional fit to the IRAS power spectrum found by Fisher et ah (1993). 
The noise term in the denominator of the Wiener filter is independent of k (cf., Fisher et al. 1993; SW, § 5.3); 
however, it is dependent on the density of galaxies, which is a decreasing function of distance in the flux- limited 
IRAS sample. As explained in Sigad et al. (1997), we therefore calculate a series of Wiener-filtered density 
fields for different noise levels, and interpolate between them to find the appropriate density field at any given 
distance. 

We wish to calculate the quadrupole due to the error in the derived density field, i.e., that due to the 
difference between Eq. ( |B10| ) and &r(k). If we expand the density field in Eq. (B4) into its Fourier components, 



substitute this difference for each component, and square the result, we find the rms contribution to uq due to 
the Wiener filter: 

UQ,Wiener) = J J r f r | P 2(»l) P M exp[i(kx • V X - k 2 • r 2 )]x 
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( [(h(k u n) - l)£r(ki) + h(h, n)e(fci)] [(/i(fc 2 , r 2 ) - l)5 T (k 2 ) + /i(A: 2 , r 2 )e(A&)] ) . (B12) 

This rather horrific expression can be simplified by multiplying out the term in brackets, realizing that the 
cross-terms vanish and that ^&r(ki)&r(k2)^ = (2 7r ) 3 -Pr(^)'5D(ki — k 2 ), where Pr(k) is the true underlying 
power spectrum, not necessarily the same as that assumed in Eq. (|B11[) . We then get two terms, one depending 
on the power spectrum, and the other due to shot noise. For the first term, the integrals over ri and r 2 separate 
to give: 

(u 2 Q,Wiener(r)} = J d 3 kP T (k) W 2 A (kR) + (u^ shot ) , (B13) 

where the new window function is given by 

W A (k,R 1 ,R) = -2 [ dr^^-[h(k,r)-l)]; (B14) 
jRi r 

compare with Eq. flB8|). We integrate from the outer volume of our peculiar velocity sample, Ri = 3000km s -1 , 
to R = 12, 800 kms -1 ; at smaller radii, the contribution to the quadrupole goes like r~ 2 , not r, and this is not 
included in our modelling of the quadrupole (Eq. |l9|). The contribution to the quadrupole from this term is 
between 1.5 and 3%, depending on which model we take for the true power spectrum. This is pleasingly close 
to the value we find for the real universe. The mock catalogs have a power spectrum set by the observed IRAS 
power spectrum (of course, with a cutoff at k < 2ir/L), and thus give a somewhat smaller contribution to this 
integral, about 1%. 

Let us now calculate the shot noise contribution to the quadrupole. It is given by: 

^Q.shot) = ^M [ «! 3 kid 3 k 2 / ^l^P^^^exp^ki-n-k,-^)] ri)e(fci)fc(fc2,^)€(fe)). 
' \ Z7T ) J J r i r 2 (B15) 

Notice now the dependence on j3, not f2; here we will make no reference to a COBE-normalized power spectrum. 
The Fourier modes are calculated in a box of side L = 25,600 kms -1 , and therefore are uncorrelated for 
Ak > 2ir/L. Thus we can write the product of the two shot noise terms as a Dirac delta function: 

2vr\ 3 „ „ . x /2tt\ 3 „ „ . x f ,o 1 



(e(k 1 )e(k 2 )) = ^(k0) ^— j 8 D 0* - k 2 ) = i^—j «5 D (k! - k 2 ) J d*r—- ; (B16) 

the expression for (e 2 (k)) comes from Fisher et al. (1993). When we insert Eq. ( B16| ) into Eq. ( |B15 ), the latter 
simplifies dramatically. The integrals over ri and r 2 now separate, giving: 

( rl )f flfUff j3 , 1 , „ . 



(«4*.W) = ^Ju{J w ' Uk) • (B17) 

where the shot noise window function looks very similar to what we have seen before: 

W shot (k,R 1 ,R) = -2 f R dr J -^^h(k,r). (B18) 
JRi r 

Notice that unlike the previous calculation, this result is independent of the true power spectrum. If we calculate 
this using the observed IRAS selection function, integrating from 3000 kms -1 to 12,800 kms -1 , we find an rms 
quadrupole of r" 1 ^ >shot (r)) 1 / 2 = 1.7(3%. 

We conclude that the 3.3% quadrupole found for the real data can be understood as a combination of the 
three effects discussed here: power on scales larger than the IRAS sample, the Wiener suppression factor, and 
shot noise; the Wiener suppression factor is the dominant one of the three. For the mock catalogs, we still do 
not completely understand why the measured residual quadrupole (< 1%) is smaller than we have calculated 
(-2%). 
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C. Properties of the Statistic x| 

In §H], we introduced the statistic xf (Eq. |2^) as a measure of the coherence of the residual field between the 
IRAS and TF data. Here we demonstrate that it has approximately the properties of a true x 2 statistic, and 
indicate how and why it departs from true x 2 behavior. 

The measure of residual coherence at separation r is 

£(r)= E <W<W ( C1 ) 

i<j 

where dy is the separation in ffiAS-distance space between objects i and j, and 5 m is the normalized magnitude 
residual, Eq. (|2~3|). The sum runs over the N p (r) distinct pairs of objects with separation r ± At; note that a 
given object may appear in more than one of these pairs. The hypothesis we wish to test is that the IRAS-TF 
residuals are incoherent, which signifies a good fit on all scales. A formal statement of this condition is that the 
individual b~ m ^ are independent random variables. Furthermore, the 5 m have been constructed to have mean 
zero and unit variance. Thus, our hypothesis of uncorrelated residuals implies that the expectation value of the 
product SmjSmj vanishes for i / j, and that the expectation value of its square is unity. 



It follows that 



E[£(t)]= E £(<WW = o- (C2) 



%<j 

=t±At 



The variance of £(t) is 

E[f(r)}= V V E{6 mii 8 m>j 8 m>k 8 mil ) . (C3) 



E 

i<j k<l 
=t±At d k i=T±Ar 



Now, the expectation value within the sum will vanish under our assumption of uncorrelated residuals unless 
i = k and j = I. (Notice that we cannot have i = I and j = k because of the ordered nature of the summation.) 



Thus, the only nonzero terms in Eq. (C3) are identical pairs, and it follows that £7[£ 2 (r)] = N p (r). 

Because £(t) is the sum of N p (r) random variables each of zero mean and unit variance, we are tempted 
to suppose that, by the central limit theorem, its distribution is Gaussian with mean zero and variance N p (r) 
when N p (t) is large. Indeed, for the 200 kms -1 bins used in its construction (cf. § |5.2] ), N p is typically ^ 10 4 . 
And, as shown in the previous paragraph, £(t) does indeed have mean zero and variance N p (t). One may also 
ask about the correlation among the £(t) for different r. Specifically, one may compute 

E[Z{n)tto)]= E E E{S m>i S m>j S m>k S m>l ) . (C4) 

i<j k<l 
dij=Ti±Ar d M =T2±Ar 

Now, it is possible to have i = k within this sum. However, because t\ 7^ T2, if i = k then j ^ I. Similarly, 
one may have j = I, but in that case i 7^ k. Thus, all of the individual expectation values in the sum vanish, 
and we find _E[£(ti)£(t2)] = 0. To the extent the above considerations hold, the £(r,) are independent Gaussian 
random variables of variance N p (ri). It then follows that the statistic x| is distributed like a x 2 variable with 
M degrees of freedom. This is the statistic proposed in the main text as a measure of goodness of fit. 

However, the central limit theorem applies only to sums of independent random variables. The individual 
products Sm^mj which enter into £(r) are uncorrelated in the specific sense E{5 m) i5 m ^)E{5 m ^5 m) i) = 8^5^ 
(where 5 K is the Kronecker-delta symbol). However, they are not strictly independent from one another. This 
is because the same object can occur in more than one pair at a given r. We thus expect the central limit to 
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apply only approximately, and the £(r) as a result are not strictly Gaussian. As a result, %| cannot be a true 
X 2 statistic. 

Furthermore, just as a single object appears in many pairs at a given r, it can appear in pairs at different 
r as well. Suppose object i contributes to both £(ti) and £(t2). Then the latter are not strictly independent, 
even though the expectation value of their product vanishes, as shown above. This factor, too, will result in a 
departure from x 2 behavior. 
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TABLE 1 



Comparison of True Parameters with Means from velmod Analyses of Mock Catalogs 



Quantity Input Value Mock Results 3, Typical Error 13 



Pi 


1.0 


0.984 ± 0.017 




0.08 


a v 


147 


149 ± 5 


20 


km s~ 1 


w LG,x° 


89 ± 8 


77 ± 12 


54 


km s~ 1 


WLG,y C 


-51 ± 10 


-50 ± 14 


63 


km s — 1 


w LG,z° 


-57 ± 9 


-55 ± 10 


45 


km s — 1 


t>A82 


10.0 


10.12 ± 0.08 




0.36 


A A82 


-13.40 d 


-13.44 ± 0.02 




0.09 


CT TF,A82 


0.45 


0.460 ± 0.006 




0.026 


b MAT 


6.71 


6.68 ± 0.05 




0.22 


A-MAT 


-5.86 d 


-5.92 ± 0.02 




0.09 


ff TF,MAT 


0.42 


0.419 ± 0.003 




0.013 



a Errors given arc in the mean. 
^Errors in a single realization. 

c Cartesian coordinates defined by Galactic coordinates. 

d Thcsc true zero points differ from those reported by Kolatt 
et ai. (1996), Table 1, because they measured distances in Mpc, 
whereas we use km s 1 . 



TABLE 2 

Numerical Results from velmod analysis of Real Data 



Quantity 


Value 


Comments 


Vq(1,1) 


37 kms -1 


at 2000 kms -1 ; cf. Eq. |l s| 


Vq(2,2) 


36 km s — 1 




Vq(1,2) 


15 kms -1 




Vq(1,3) 


113 kms -1 




Vq(2,3) 


-24 kins" 1 




u v 


125 kms" 1 




w LG,x 


-30 km! -1 




W LG ,y 


-10 kms" 1 




w LG,z 


30 kms -1 




b A82 


10.36 ± 0.36 


10.29 ± 0.22 (Mark III value) 


A A82 


-5.96 ± 0.09 


-5.95 ± 0.04 (Mark III value) 


°"TF,A82 


0.464 ± 0.026 


0.47 ± 0.03 (Mark III value) 


b MAT 


7.12 ± 0.22 


6.80 ± 0.08 (Mark III value) 


A MAT 


-5.75 ± 0.09 


-5.79 ± 0.03 (Mark III value) 


ff TF,MAT 


0.453 ± 0.013 


0.43 ± 0.02 (Mark III value) 


01 


0.492 ± 0.068 


With Quadrupolc 


01 


0.563 ± 0.074 


Without Quadrupolc 


01 


0.489 ± 0.084 


A82 data only 


01 


0.498 ± 0.107 


MAT data only 


01 


0.453 ± 0.093 


< cz LG < 1350 kms -1 


01 


0.495 ± 0.133 


1350 < cz LG < 2150 kms -1 


01 


0.573 ± 0.142 


2150 < cz L G < 3000 kms -1 


01 


0.521 ± 0.050 


wlg — 0; &v fixed to 250 kms -1 


01 


0.491 ± 0.045 


wlg — 0; &v fixed to 150 kms - 


01 


0.544 ± 0.071 


With Quadrupolc; 500 kms -1 smoothing 


01 


0.635 ± 0.083 


Without Quadrupolc; 500 kms" 1 smoothing 


01 


0.510 ± 0.038 


TF parameters fixed at Mark III values; with quadrupolc 


01 


0.517 ± 0.039 


TF parameters fixed at Mark III values; without quadrupolc 



We did not do a likelihood search in parameter space to find formal error bars on quantities 
other than /3j . Error estimates fopi the TF parameters come from averaging over the mock 
catalog VELMOD runs; sec Tabic H. 
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